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Introduction: 

A  leading  cause  of  cancer  deaths  among  women  worldwide  is  due  to  breast 
cancer  and  new  therapies  and  optimal  treatments  are  continuously  being  conceived  and 
explored  to  better  control  or  even  cure  this  disease.  (1).  Diagnosis  and  therapeutic 
management  of  the  breast  tumor  remain  significant  medical  challenges,  hence  early 
detection,  diagnosis,  and  timely  treatments  are  essential  to  successful  health  care  (2). 
Magnetic  resonance  imaging  and  spectroscopy  (MRI/MRS)  have  gained  significant 
importance  during  the  last  fifteen  years  for  the  diagnosis  and  monitoring  of  breast  cancer 
therapy.  The  sensitivity  of  MRI/MRS  for  anatomical  delineation  is  very  high  and  the 
consensus  is  that  MRI  is  more  sensitive  in  detection  than  x-ray  mammography. 
Advantages  of  MRS  include  delivery  of  biochemical  information  about  tumor 
metabolism,  which  can  potentially  assist  in  the  staging  of  cancers  and  monitoring 
responses  to  treatment.  High  sensitivity  is  a  major  advantage  of  contrast  enhanced  MRI, 
but  its  diagnostic  relevance  in  the  future  will  largely  depend  on  improvements  in 
specificity.  Current  approaches  in  the  application  of  MRI  to  breast  tumors  aim  to 
improve  specificity  and  sensitivity  (4-17).  Increased  specificity  is  necessary  to  reduce  the 
number  of  biopsies  performed  to  confirm  false  positive  findings.  Diffusion-weighted 
imaging  (DWI)  is  another  MR  based  technique  that  probes  the  microstructure  of  tissues 
and  is  sensitive  to  the  degree  to  which  motion  of  water  molecules  is  restricted  in  relation 
to  how  packed  together  cells  are  (17,  18).  It  has  been  reported  that  high  resolution  DWI 
may  add  valuable  functional  information  to  conventional  MR  protocols  with  short 
measurement  times  for  the  diagnosis  of  breast  cancer  and  improve  the  specificity  of  MR 
imaging  (19-21).  However,  new  technological  developments  are  necessary  to  assess  their 
role  in  breast  diagnosis.  A  method  capable  of  identifying  biochemical  characteristics  non- 
invasively  in  the  tumor  lesions  that  can  be  used  in  conjunction  with  MRI  is  proton  ('H) 
MR  Spectroscopy  (MRS).  Researchers  have  shown  that  1 H  MRS  can  be  used  to 
characterize  breast  cancers  with  improved  diagnostic  accuracy  (22-26).  Unfortunately, 
multi-voxel  based  novel  MR  spectroscopic  imaging  (MRSI)  techniques  using  the  speed 
advantage  offered  by  echo-planar  imaging  (EPI)  and  improved  spectral  resolution  offered 
by  two-dimensional  (2D)  MR  spectroscopy  (MRS)  have  not  been  fully  explored  in  breast 
cancer  studies  so  far.  Hence,  a  major  task  of  this  project  is  to  combine  echo-planar 
correlated  spectroscopic  imaging  (EP-COSI)  with  DWI  approach  for  improving  the 
overall  specificity  of  breast  cancer  detection. 

Body: 

i)  Proposed  Task  1:  To  further  optimize  the  multi-voxel  based  extension  of  the  correlated 
spectroscopy  (COSY)  sequence,  in  which  two  spectral  encodings  will  be  combined  with 
two  spatial  encodings.  This  four-dimensional  (4D)  data  acquisition  scheme  will  be 
accomplished  utilizing  the  echo-planar  imaging  (EPI)  approach  that  is  commonly  used 
for  spatial  encoding  in  MRI  including  DWI.  (Months  1-6). 

This  task  was  completed  already  (September  2010-May  201 1)  as  reported  previously. 

ii)  Proposed  Task  2:  To  evaluate  the  EP-COSI  data  using  a  breast  phantom  containing 
two  concentric  spheres,  the  inner  one  containing  several  metabolites  which  have  been 
reported  in  breast  tissues  surrounded  by  the  outer  phantom  containing  corn  oil  to  mimic 
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fatty  tissues  known  to  be  in  breast  tissues,  and  to  optimize  the  echo  speed-factor  and 
other  acquisition  parameters  using  the  phantom  (Months  1-6). 

This  task  was  completed  already  (September  2010-August  201 1)  as  reported  in  our  2012 
report. 


nil  Proposed  Task  3:  To  develop,  evaluate  and  optimize  the  prior-knowledge  basis  set 
spectra  using  the  GAMMA-simuIation  and  breast  phantom  solutions  as  prior  knowledge 
for  the  multi-voxel  based  COSY  spectra  recorded  using  the  3T  MRI  scanner  (Months  3- 

9f 

Using  the  GAMMA  library  (27),  a  prior-knowledge  basis-sets  containing  metabolites  and 
lipids  have  been  prepared  as  reported  previously. 

iv)  Proposed  Task  4:  To  record  the  EP-COSI  spectra  in  the  fatty,  glandular  and  ductal 
areas  of  healthy  breasts.  Twenty  healthy  female  volunteers  (25-70  years  old)  with  no 
previous  history  of  breast  cancer  will  be  investigated.  (Months  6-24). 

Eleven  healthy  women  (age=26-59  years)  without  any  prehistory  of  breast  cancer  were 
investigated  using  the  DWI  and  EP-COSI  MRI  protocol  during  the  first  two  years  as 
reported  previously.  Six  more  healthy  women  were  investigated  during  the  current  year. 

v)  Proposed  Task  5:  To  record  the  multi-voxel-based  2D  spectra  in  patients  with  benign 
and  malignant  breast  cancer.  The  breast  metabolite  and  lipid  concentrations  calculated 
from  the  multi-voxel  data  using  the  ProFit  algorithm  will  be  compared  with  LC-Model 
processed  ID  spectral  based  MRS1  data.  Twenty  patients  with  biopsy-proven  breast 
cancer  (ductal  carcinoma  and  invasive  lobular  cancer),  twenty  patients  with  benign 
breast  tumor  (fibroadenoma,  proliferative  fibrocystic  change  and  papillomas)  will  be 
investigated  (Months  6-24). 

We  had  reported  studying  three  women  with  benign  (age=28-44  years)  and  2  with 
malignant  breast  cancer  (age=6 1  years)  in  our  last  annual  report.  In  the  current  year, 
using  a  dedicated  breast  MRI  coil,  we  have  investigated  four  malignant  and  three  benign 
cancer  patients  using  the  4D  EP-COSI  sequence  on  the  3T  MRI  scanner. 

Shown  in  Fig.  1  are  the  following:  A)  An  axial  MRI  slice  image  recorded  in  a 
38yo  malignant  breast  cancer  patient  showing  the  EP-COSI  localization  (white  box).  B) 
an  extracted  COSY  spectrum(2ml)  recorded  from  the  affected  breast  lesion  and  C) 
another  extracted  COSY  spectrum(2ml)  recorded  from  the  unaffected  breast  fatty  region. 
As  reported  previously  (28),  the  2D  COSY  spectrum  from  the  suspicious  mass  showed 
elevated  water  and  decreased  lipid  resonances. 

Fig.  2  shows  the  following:  A)  An  axial  MRI  slice  image  recorded  in  a  3  lyo 
benign  breast  cancer  patient  showing  the  EP-COSI  localization.  B)  Multi-voxel  EP-COSI 
spectra;  C)  and  D)  show  2ml  COSY  spectra  recorded  from  the  affected  mass  and  a 
healthy  region.  These  results  are  also  in  agreement  with  our  earlier  data  on  the  SV-based 
L-COSY  sequence  (29-30). 
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Figure  1.  (A)  An  axial  MRI  slice  image  recorded  in  a  38  yo  subject  with  malignant 
breast  cancer  showing  the  EP-COSI  localization;  the  small  white  box  shows  the  regions 
of  interest  (ROI)  localized  the  EP-COSI  sequence  and  the  large  white  box  with  16x16 
grids  shows  the  spatial  localization  encoded  by  the  EP-COSI  sequence.  (B)  A  2ml  2D 
COSY  spectrum  from  the  affected  breast  region.  (C)  A  2ml  spectrum  showing  healthy 
fatty  breast  region. 
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Figure  2.  (A)  An  axial  MRI  slice  image  recorded  in  a  3 1  yo  subject  with  benign  breast 
cancer  (a  hypo-intense  mass)  showing  the  EP-COSI  localization;  the  small  white  box 
shows  the  regions  of  interest  (ROI)  localized  the  EP-COSI  sequence  and  the  large  white 
box  with  16x16  grids  shows  the  spatial  localization  encoded  by  the  EP-COSI  sequence. 
(B)  An  extracted  2D  COSY  spectrum  from  the  affected  breast  regions.  (C)  2ml  spectrum 
from  the  fatty  breast  region. 


7 


vi)  Proposed  Task  6:  To  record  multi-slice  DWI  in  twenty  patients  with  biopsy-proven 
breast  cancer,  twenty  patients  with  benign  breast  tumor  and  twenty  healthy  women  and  to 
calculate  the  ADC  maps.  (Months  6-24). 

Figure  3.  An  axial  ADC  slice  image  recorded  in  the  38  yo  subject  with  a  malignant  mass 
using  the  DWI  sequence. 


Figure  4.  An  axial  ADC  slice  image  recorded  in  the  31  yo  subject  with  a  benign  mass  in 
both  breasts  using  the  DWI  sequence. 


V  I 
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The  diffusion  weighted  imaging  (DWI)  has  been  recorded  in  three  women  with  benign 
and  4  with  malignant  breast  cancer  3T  MRI  scanner  using  a  dedicated  breast  MRI  coil. 

vii )  Proposed  Task  7:  To  correlate  the  EP-COSI  findings  with  that  of  DWI  in 
differentiating  benign  from  malignant  breast  cancers,  and  to  calculate  specificity, 
sensitivity  and  accuracy  of  the  MRSI  and  DWI  data  in  differentiating  benign  from 
malignant  tumors.  (Months  6-24). 
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This  task  will  be  completed  during  the  no-cost  extension  period. 

Key  Research  Accomplishments 

•  The  4D  EP-COSI  has  been  successfully  evaluated  in  17  healthy  women,  5  benign 
and  6  malignant  breast  cancer  patients  so  far  using  the  UCLA  Radiology  Siemens 
3T  MRI  scanner  equipped  with  a  dedicated  breast  phased-array  assembly.  As 
summarized  in  the  previous  year  reports,  the  4D  EP-COSI  sequence  is  available  at 
UCLA  only  as  the  sequence  is  not  supplied  by  any  of  MRI  manufacturers. 

•  The  DWI-MRI  protocol  was  successfully  evaluated  in  a  total  of  17  healthy 
women,  5  benign  and  6  malignant  breast  cancer  patients  so  far. 

•  Our  pilot  results  so  far  clearly  demonstrate  that  the  EP-COSI  spectroscopic 
imaging  sequence  can  be  combined  with  a  clinical  breast  DWI  protocol  with  the 
total  duration  of  less  than  an  hour.  The  protocol  is  completely  safe  enough  to  be 
included  in  any  MRI  protocol  to  be  evaluated  in  breast  cancer  for  improving  the 
overall  specificity. 

•  Our  recent  work  has  focused  on  testing  retrospective  Maximum  Entropy  and 
Compressed  Sensing  of  the  4D  EP-COSI  data  so  that  the  acquisition  can  be 
accelerated  to  less  than  10  minutes  which  will  lead  to  less  imaging  time  for 
patients  with  breast  cancer. 


Reportable  Outcomes: 

A.  Peer-reviewed  Publications:  1)  Bums  B,  Wilson  N,  Furuyama  JK  and  Thomas  MA. 
Non-uniformly  under-sampled  multi-dimensional  spectroscopic  imaging  in  vivo: 
maximum  entropy  versus  compressed  sensing  reconstruction.  NMR  Biomed  20 14;27 : 
191-201. 

2)  Burns  B,  Wilson  N  and  Thomas  MA.  Split-Bregman-based  Group-Sparse 
Reconstruction  of  multi-dimensional  MR  Spectroscopic  Imaging  data.  ISBI  Beijing,  2014 
(in  press) 

B.  Presentations:  We  report  one  conference  presentation  as  shown  below. 

•  An  oral  presentation  at  the  ISBI  conference,  Beijing.  China,  April28-May  2,  2014  (to 
be  presented) 


C.  Books:  None  on  Breast  Cancer  Research  based. 

Conclusions:  So  far,  the  MR  scanning  protocol  including  EP-COSI  and  DWI-MRI 
has  been  evaluated  in  17  healthy  women,  5  benign  and  6  malignant  breast  cancer.  We 
will  continue  to  recruit  8-10  malignant,  4-6  benign  and  3  healthy  women  during  the  next 
year.  The  UCLA  mammo  co-investigators,  namely  Dr.  Naneltte  DeBruhl  and  Dr. 
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Lawrence  Bassett  have  promised  us  to  help  recruiting  more  patients  in  the  current  no-cost 
extension  year. 
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Non-uniformly  under-sampled  multi-dimensional 
spectroscopic  imaging  in  vivo :  maximum  entropy 
versus  compressed  sensing  reconstruction 


Brian  Burnsa,b/ 
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The  four-dimensional  (4D)  echo-planar  correlated  spectroscopic  imaging  (EP-COSI)  sequence  allows  for  the  simultaneous 
acquisition  of  two  spatial  (ky,  kx)  and  two  spectral  (f2,  t-,)  dimensions  in  vivo  in  a  single  recording.  However,  its  scan  time  is 
directly  proportional  to  the  number  of  increments  in  the  ky  and  fn  dimensions,  and  a  single  scan  can  take  20-40  min 
using  typical  parameters,  which  is  too  long  to  be  used  for  a  routine  clinical  protocol.  The  present  work  describes  efforts 
to  accelerate  EP-COSI  data  acquisition  by  application  of  non-uniform  under-sampling  (NUS)  to  the  ky-t-,  plane  of 
simulated  and  in  vivo  EP-COSI  datasets  then  reconstructing  missing  samples  using  maximum  entropy  (MaxEnt)  and 
compressed  sensing  (CS).  Both  reconstruction  problems  were  solved  using  the  Cambridge  algorithm,  which  offers 
many  workflow  improvements  over  other  /,-norm  solvers.  Reconstructions  of  retrospectively  under-sampled 
simulated  data  demonstrate  that  the  MaxEnt  and  CS  reconstructions  successfully  restore  data  fidelity  at  signal-to-noise 
ratios  (SNRs)  from  4  to  20  and  5x  to  1.25x  NUS.  Retrospectively  and  prospectively  4x  under-sampled  4D  EP-COSI 
in  vivo  datasets  show  that  both  reconstruction  methods  successfully  remove  NUS  artifacts;  however,  MaxEnt  provides 
reconstructions  equal  to  or  better  than  CS.  Our  results  show  that  NUS  combined  with  iterative  reconstruction  can  reduce 
4D  EP-COSI  scan  times  by  75%  to  a  clinically  viable  5  min  in  vivo,  with  MaxEnt  being  the  preferred  method.  Copyright  © 
2013  John  Wiley  &  Sons,  Ltd. 

Keywords:  EP-COSI;  maximum  entropy;  compressed  sensing;  non-uniform  under-sampling;  spectroscopy;  spectroscopic 
imaging 


INTRODUCTION 

Changes  in  metabolite  concentrations  as  a  result  of  the  altered 
metabolism  of  cancer  can  be  detected  non-invasively  using 
one-dimensional  (ID)  MRS  in  vivo  (1-3).  However,  the  overlap 
of  spectral  peaks  in  ID  MRS  is  a  major  impediment  to  the 
identification  of  individual  metabolites.  Two-dimensional  (2D) 
MRS  has  increased  spectral  dispersion  over  ID  MRS  and  can 
disentangle  overlapping  complex  spectral  peaks  (4).  Single-voxel 
2D  MRS  has  been  shown  to  increase  the  specificity  and  sensi¬ 
tivity  of  tumor  grade  classification  when  used  with  dynamic 
contrast-enhanced  MRI  in  the  breast  (5).  However,  the  acquisi¬ 
tion  of  multiple  f-i  increments  per  voxel  to  form  the  second 
spectral  dimension  limits  its  ability  to  provide  multi-voxel 
coverage  because  of  the  long  scan  times  needed  to  combine 
two  spectral  and  two  spatial  dimensions. 

With  the  advent  of  echo-planar  spectroscopic  imaging  (EPSI), 
MRSI  scans  with  one  spectral  and  two  spatial  dimensions  can 
be  completed  within  clinically  acceptable  times  by  interleaving 
the  acquisition  of  a  spatial  and  spectral  dimension  within  the 
EPSI  readout  (6-8).  The  four-dimensional  (4D)  echo-planar  corre¬ 
lated  spectroscopic  imaging  (EP-COSI)  (9)  sequence  allows  the 
acquisition  of  two  spatial  (ky,  kx)  and  two  spectral  (t2,  fq)  dimen¬ 
sions  in  a  single  recording  to  form  4D  MRSI.  The  sequence  inter¬ 
leaves  the  acquisition  of  the  kx  and  t2  dimensions  within  the  EPSI 
readout,  but  ky  and  f-,  are  incrementally  acquired  as  indirect 
dimensions  during  each  TR.  The  EP-COSI  sequence  has  the  ben¬ 
efits  of  increased  spectral  dispersion  and  multi-voxel  support, 
which  improves  metabolite  identification  over  multiple  spatial 


regions  simultaneously;  however,  its  scan  time  is  directly  propor¬ 
tional  to  the  number  of  increments  in  the  ky  and  f,  dimensions. 
An  EP-COSI  scan  using  typical  parameters  of  TR/TE  =  1 .5  s/30  ms 
and  /cy/fq  =  16/100  can  take  40  min,  which  is  too  long  to  be  used 
within  a  routine  clinical  protocol. 
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When  the  Fast  Fourier  Transform  (FFT)  is  used  to  transform 
uniformly  sampled  4D  MRSI  data  ( ky ,  kx,  t2,  fq)  to  the  spatial, 
spectral  domain  (V,  X,  F2,  Fq),  decreasing  scan  times  require  a 
reduction  in  either  the  ky  spatial  or  fq  spectral  dimension  through 
truncation  or  lower  sampling  rates,  and  a  corresponding  unwanted 
reduction  in  resolution  or  bandwidth.  Flowever,  non-uniform 
under-sampling  (NUS)  of  the  spatial,  spectral  /(y-fq  plane,  in  combi¬ 
nation  with  iterative  non-linear  reconstruction,  can  be  used  to 
accelerate  the  collection  of  4D  MRSI  data  in  vivo,  whilst  preserving 
the  spatial  and  spectral  resolutions  and  bandwidths  (10). 

Earlier  work  has  demonstrated  the  feasibility  of  under-sampling 
the  mixed-domain  /cy-fq  plane  of  a  4D  echo-planar  ./-resolved 
spectroscopic  imaging  (EP-JRESI)  dataset  and  reconstructing  the 
missing  points  with  compressed  sensing  (CS)  (10),  a  popular 
method  of  non-linear  iterative  image  reconstruction  which 
promotes  data  sparsity  in  the  reconstruction  domain  and  data 
fidelity  in  the  sample  domain  (1 1,1 2).  The  nature  of  spatial,  spectral 
NUS  artifacts  in  the  /(y-fq  plane  was  explored,  and  it  was  shown 
that/, -norm-based  CS  reconstruction  is  a  viable  means  of  reducing 
the  scan  times  of  4D  EP-JRESI  in  vivo  through  NUS.  In  recent  years, 
CS  reconstruction  has  been  successfully  applied  to  NUS  MRI 
(13,14),  three-dimensional  (3D)  MRSI  (15),  dynamic  MRI  (16,17), 
and  multi-dimensional  Nuclear  Magnetic  Resonance  (NMR)  (23). 

Maximum  entropy  (MaxEnt)  image  reconstruction  is  an  alterna¬ 
tive  non-linear  iterative  reconstruction  technique  to  CS.  Rather 
than  minimizing  transform  sparsity,  it  maximizes  the  entropy  of 
the  data  in  the  reconstruction  domain,  whilst  preserving  data 
fidelity  in  the  sample  domain  (1 8,1 9).  MaxEnt  has  been  successfully 
used  to  reconstruct  under-sampled  images  in  astronomy  and 
multi-dimensional  spectra  in  NMR  (19-21).  Flowever,  MaxEnt  has 
not  been  applied  to  the  mixed-domain  /cy-fq  plane  of  a  4D  MRSI 
dataset  in  vivo. 

The  use  of  entropy  as  a  regularizer  in  image  reconstruction 
predates  /,-norm-based  CS  reconstruction  and  continues  to  be 
used  extensively  in  the  reconstruction  of  under-sampled  NMR 
spectra  in  spite  of  the  popularity  of  CS  in  other  fields.  It  was  first 
suggested  by  Frieden  (18)  in  the  early  1970s  after  Jaynes  (22) 
proposed  the  idea  of  the  Principle  of  Maximum  Entropy,  which 
describes  the  MaxEnt  distribution  as  the  'maximally  non-committal 
distribution  with  regard  to  unavailable  data'.  This  principle  pre¬ 
sents  the  MaxEnt  prior  as  one  that  assumes  nothing  about  the 
unavailable  data;  by  assuming  nothing  about  those  points,  their 
possible  values  are  all  equally  likely  to  occur,  and  the  reconstruc¬ 
tion  is  that  which  most  closely  conforms  to  the  uniform  distribu¬ 
tion,  i.e.  is  flat.  Peaks  in  the  reconstruction  domain  are  the  result 
of  signals  from  the  sampled  data,  and  any  artifacts  from  the  miss¬ 
ing  data  points  are  removed  because  they  represent  states  of  low 
entropy  that  are  not  the  result  of  k- space  or  time-domain  signals. 

In  this  article,  we  compare  the  MaxEnt  and  /, -norm-based  CS 
reconstructions  of  NUS  4D  EP-COSI  data  and  show  that 
MaxEnt  is  a  viable  alternative  to  CS  for  reducing  scan  times  4x 
in  human  breast  in  vivo.  We  quantitatively  characterize  the 
MaxEnt  and  CS  reconstructions  by  comparing  results  for 
retrospectively  NUS-simulated  4D  EP-COSI  data  at  varying  levels 
of  signal-to-noise  ratio  (SNR)  and  NUS  rates.  We  show  that  retro¬ 
spectively  4x  NUS  4D  EP-COSI  in  vivo  breast  data  reconstructed 
using  either  MaxEnt  or  CS  show  a  comparable  spatial,  spectral 
resolution  to  the  fully  sampled  data.  In  addition,  we  show  that 
MaxEnt  and  CS  reconstructions  of  prospective  4x  NUS  4D 
EP-COSI  scans  from  the  same  breast  study  as  the  retrospective 
data,  using  the  same  mask  and  sequence  parameters,  compare 
favorably  with  the  retrospective  and  fully  sampled  data. 


Throughout  this  article,  the  NUS  dataset  that  has  not  been 
reconstructed  and  has  zeros  in  place  of  missing  samples  is 
referred  to  as  the  zero-augmented  dataset  to  distinguish  it  from 
the  MaxEnt  and  CS  reconstructions. 

EXPERIMENTAL  DETAILS 

4D  MaxEnt  and  CS  reconstruction:  theory 

MaxEnt  and  /, -norm-based  CS  were  used  to  reconstruct  the  NUS 
ky-t-]  plane  of  the  4D  EP-COSI  datasets. 

/, -norm-based  CS  image  reconstruction  of  4D  MRSI  data  is 
formulated  as  a  constrained  convex  optimization  problem 
(10-12): 

minimize  ||  (/rm  ||  q  s.t  \\KTm  —  d\\2<C0  [1] 

where  m  =  (y,x,  F2,  Fq)  is  the  reconstructed  spatial,  spectral-do- 
main  data,  T  is  the  4D  Fourier  operator,  K  is  the  NUS  mask  that 
determines  which  samples  were  acquired  in  the  /cy-fq  plane, 
d  =  (ky,k„t2,ti)  is  the  /c-space,  time-domain  sampled  data,  C0  is 
the  standard  deviation  of  the  noise  in  d  and  y/  is  a  known  sparse 
transform,  y/  was  chosen  to  be  the  identity  transform  because  m 
was  already  self  sparse  as  shown  in  ref.  (10). 

MaxEnt  image  reconstruction  of  4D  MRSI  data  solves  a  similar 
problem  to  CS,  but  uses  Sq/2  entropy  instead  of  the  lrnorm  (19-21): 

maximize  S(m)-,/2  s.t  \\KTm  -  d\\j<C0  [2] 

where  5(m), ,  is  the  entropy  of  the  estimated  spectrum,  and  the 
remaining  terms  are  identical  to  those  in  the  CS  problem.  Sq/2 
entropy  is  a  concave  function  with  a  global  maximum  and  no  local 
extrema,  and  so  there  is  a  single  solution  that  satisfies  the  problem 
within  the  feasible  set  of  solutions  defined  by  the  data  fidelity 
constraint  as  shown  in  Fig.  1  (24).  As  can  be  seen,  S1/2  entropy 
has  a  global  extremum  and  has  slightly  more  curvature  than  the 
lrnorm,  but  has  far  less  curvature  than  the  l2-norm. 

The  entropy  used  in  the  MaxEnt  reconstruction  was  not  the 
often  used  —  Iplog(p)  entropy  introduced  by  Shannon  (25), 
but  the  Sq/2  entropy  derived  by  Daniell  and  Flore  (20)  specifically 


Figure  1.  Plot  showing  -S1/2  entropy  (blue),  lrnorm  (red)  and  l2-norm 
(green).  Each  function  has  been  normalized  to  equal  unity  at  |x|  =  1. 
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for  NMR  spectra  originating  from  spin 1  /2  nuclei,  such  as  ’H,  used 
in  MRSI: 


diagonal  blocks  of  the  ||m||n  and  S(m,)i  Hessians  were  stored 
in  memory  during  reconstruction,  not  the  full  matrices. 


where  def  is  a  scaling  parameter  related  to  the  sensitivity  of  the 
scanner  and  is  calculated  for  m  of  length  N  as  \/Wn  (21).  S1/2  is 
used  because  the  underlying  physical  processes  that  produce  an 
MR  spectrum  are  not  based  on  discrete  particle  events  and  so 
cannot  be  modeled  by  simple  Poisson-distributed  processes  as 
required  for  the  derivation  of  Shannon  entropy  (26).  They  are 
governed  by  the  density  matrix  of  the  spin  system  under  inves¬ 
tigation,  and  so  the  statistical  distribution  is  different.  Equation 
[3]  was  derived  from  first  principles  using  both  a  classical  spin 
model  and  a  quantum  mechanical  model.  Neither  model  made 
any  assumptions  on  the  initial  state  of  the  spin  system  nor  the 
pulse  sequence  used.  This  equation  can  be  applied  to  any  MR 
spectrum  originating  from  spin  1  /2  nuclei,  and  addresses  previ¬ 
ous  concerns  regarding  the  use  of  entropy  in  MRS  and  MRI 
reconstruction  (27). 

In  order  to  remove  any  differences  between  the  MaxEnt  and 
CS  reconstructions  caused  by  differences  in  the  solvers  used, 
both  problems  were  solved  by  a  Matlab  implementation  of  the 
Cambridge  algorithm  (19).  This  recasts  the  image  reconstruction 
problem  into  an  unconstrained  convex  optimization  problem 
and  uses  a  variant  of  the  conjugate  gradient  method  to  itera¬ 
tively  find  the  extrema  in  two  phases;  the  first  phase  minimizes 
the  fidelity  constraint  and  the  second  phase  minimizes  or  maxi¬ 
mizes  the  objective  function,  while  keeping  the  fidelity 
constraint  minimized.  The  stopping  criterion  for  the  problem  is 
reached  when  the  gradients  of  the  objective  0{m)  and  the 

fidelity  constraint  C(m)  are  parallel:  1^^  —  |^|  <  0.001.  Specific 

details  on  the  algorithm  and  modifications  to  accommodate 
multi-dimensional  MR  data  can  be  found  in  ref.  (21). 

The  Cambridge  algorithm  calculates  the  gradient  V||m||1eCN 
and  Hessian  V2||m||1eC2A,)<2N  of  the  objective  function,  which 
are  not  defined  for  || m  ||  n  when  m,=  0.  Therefore,  in  order  to 
solve  the  lrnorm- based  CS  reconstruction  problem,  ||m||,  was 
redefined  as: 


II  mil,  =  llJn(mi)2  +  Z(mi)2  +  e  [4] 

/=  i  ’ 

where  7?.  and  3  were  the  real  and  imaginary  components  of  m„ 
and  e  is  a  small  non-zero  value  to  prevent  m,  =  0.  The  gradient 
and  Hessian  were  then  defined  as: 

V||m||1  =  Mr1m  [5] 


V2 1|  m  ||  t 


M  1  (1  -  H(mi))\m,\  2 


-3(m,-)72(m,)|m,j  3 


.  -3(m,)72(m,)|m,|  3 


H-1(1 


3(m,-))|m,-| 


[6] 


=  2x2  block  diagonal  matrix 

where  M/eCWxW  is  a  diagonal  matrix  with  wn=  |m,|,  and  element  m , 
is  associated  with  the  /'th  2x2  block  in  the  Hessian.  The  gradient 
and  Hessian  of  S(m,) ,  are  defined  in  ref.  (21).  Only  the  2x2 


Sample  mask  generation 

The  ky-ti  plane  of  the  4D  EP-COSI  datasets  used  in  these 
experiments  was  under-sampled  using  2D  Poisson-gap  sample 
masks  that  were  generated  using  a  modified  ID  Poisson-gap 
process  (28).  In  1 D  Poisson-gap  sample  masks,  the  gaps  between 
samples  follow  a  Poisson  distribution,  whereas  the  2D  extension 
follows  the  convention  that  gaps  between  spaces  follow  a 
Poisson  distribution.  However,  both  conventions  result  in  the 
spaces  and  sample  points  following  Poisson  distributions. 
Poisson-distributed  masks  avoid  large  gaps  between  samples, 
which  are  detrimental  to  the  reconstruction,  while  ensuring  that 
the  samples  are  randomly  distributed  (29).  Compared  with  other 
distributions,  Poisson-distributed  masks  create  the  fewest 
aliasing  artifacts  in  the  Fourier  domain  and  preserve  the  SNR  of 
the  under-sampled  data  (30). 

The  effects  of  the  sample  mask  on  the  peak  amplitude  and 
lineshape  of  spectral  reconstructions  are  well  documented 
(29,31,32).  Sample  mask  densities  that  follow  the  time-domain 
NMR  signal  envelope  and  sample  more  points  at  higher  SNR 
have  spectral  reconstructions  with  lower  root-mean-square 
errors  (RMSEs)  and  less  non-linearity  compared  with  sample 
masks  that  do  not.  The  f,  dimensions  of  the  EP-COSI  datasets 
in  this  work  were  apodized  with  a  sine-squared  filter  to  enhance 
the  cross-peaks  (9),  but,  because  of  T2*  decay,  the  filtered  EP- 
COSI  data  had  a  skewed  sine-squared  signal  envelope;  therefore, 
the  2D  Poisson-gap  sample  mask  density  was  modulated  along 
t,  with  a  skewed  sine-squared  function  (33).  The  ky  dimension 
was  modulated  by  an  exponential  decay  function  similar  to  that 
used  previously  to  maximize  spatial  SNR  (34). 

The  sample  density  of  a  Poisson-gap  mask  can  be  modulated 
by  the  rate  parameter  2  which  determines  both  the  mean  and 
variance  of  a  Poisson  distribution.  The  probability  of  generating 
a  gap  g  from  a  Poisson  distribution  is  characterized  by: 


p(g,2)  =  (2s -e  l)/g!  [7] 


For  large  2,  large  values  of  g  are  more  likely,  and  for  small  2, 
small  values  of  g  are  more  likely.  Therefore,  the  probability  of  g 
can  be  modulated  by  varying  the  value  of  2  according  to  a  sine 
or  exponential  decay  function,  and  the  probability  of  large  gaps 
between  spaces  can  be  increased  where  the  SNR  is  highest  in 
the  MR  signal  envelope  (28).  To  generate  g  as  a  function  of  2,  a 
Poisson  process  can  be  simulated  using  various  techniques  that 
do  not  depend  on  an  a  priori  knowledge  of  g  as  above  (35).  For 
these  experiments,  the  poissrnd(X)  function  in  Matlab  was  used 
to  generate  g  as  a  function  of  2.  It  takes  as  input  an  array  of  2 
and  returns  an  array  of  gaps  with  local  mean  and  variance  2. 

2D  Poisson-gap  sample  masks  were  iteratively  generated  in 
Matlab  by  combining  the  ID  distributions  of  f,  and  ky  until  the 
desired  NUS  rate  was  reached. 

Examples  of  the  2D  2,  gap  and  mask  arrays  generated  by  2D 
Poisson-gap  are  shown  in  Fig.  2;  as  the  size  of  2  and  the  spacing 
gaps  increase,  the  sample  density  increases  in  that  area  of  the 
mask.  The  magnitude  point  spread  function  (PSF)  of  the  mask 
is  shown  and  demonstrates  the  viability  of  this  approach;  the 
single  dominant  central  peak  with  small  side-lobes,  surrounded 
by  low-amplitude,  incoherent  artifacts,  is  the  desired  profile  of 
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Figure  2.  Poisson  sample  mask  creation  along  the  ky-fq  plane  for  echo- 
planar  correlated  spectroscopic  imaging  (EP-COSI).  Top:  modulated 
values  of  A.  Upper  middle:  Poisson-distributed  values  for  each  A  indicat¬ 
ing  the  gap  between  spaces  in  the  sample  distribution.  Lower  middle: 
resulting  two-dimensional  (2D)  Poisson-gap  sample  mask.  Bottom:  mag¬ 
nitude  point  spread  function  (PSF)  of  the  2D  sample  mask. 


a  NUS  mask  PSF  (36,37).  Incoherent  sampling  artifacts  will  have 
low  amplitudes  and  spurious  peaks  caused  by  coherent  aliasing 
will  be  negligible. 


They  were  then  copied  into  an  8x8  spatial  grid  to  simulate 
spatially  distributed  metabolites  as  follows:  the  upper  left  quadrant 
contained  2  x2  voxels  of  tCho,  the  upper  right  quadrant  contained 
2x2  voxels  of  Cr,  the  lower  left  quadrant  contained  2x2  voxels  of 
Glx  and  the  lower  right  quadrant  contained  2x2  voxels  of  Asp. 

As  a  result  of  under-sampling  the  /cy-fq  plane,  the  spatial, 
spectral  artifacts  caused  the  tCho  and  Glx  voxels  to  alias  into 
each  other  and  the  Cr  and  Asp  voxels  to  alias  into  each  other, 
as  illustrated  by  the  zero-augmented  data  at  the  bottom  of  Fig.  3. 
This  changed  the  integrated  peak  area  contained  within  each 
metabolite  as  the  spatial,  spectral  separation  between  the 
metabolites  broke  down. 

Noise  was  added  to  the  simulated  noise-free  4D  EP-COSI 
dataset  to  model  the  SNRs  of  2-20  in  increments  of  2.  It  was 
under-sampled  5x,  2.5x,  1.67X  and  1.25x,  and  then  separately 
reconstructed  by  MaxEnt  and  CS.  The  SNR  was  varied  by 
simulating  different  levels  of  thermal  noise  in  the  dataset  by 
adding  univariate  Gaussian  noise  to  the  noise-free  real  and 
imaginary  channels  of  the  4D  EP-COSI  dataset  (40).  The  desired 
SNR  was  achieved  by  ensuring  that  the  additive  noise  signal 
power  (cr2)  was  equal  to  1/SNR  of  the  noise-free  dataset  signal 
power  (co2),  such  that: 

noisy  data  =  noise  free  data  +  - x  noise  [8] 

a2  x  SNR 


MR  simulations 

The  effects  of  SNR  and  under-sampling  rate  on  the  reconstructions 
were  quantitatively  assessed  using  a  noise-free  simulated  4D  EP- 
COSI  dataset  that  contained  choline +  glycerophosphocholine  + 
phosphocholine  (total  choline,  tCho),  glutamate  +  glutamine 
(Glx),  creatine  (Cr),  aspartate  (Asp)  or  nothing  in  each  voxel,  as 
represented  in  the  top  of  Fig.  3  by  a  diagonal  peak  from  each 
metabolite.  Each  metabolite  was  simulated  using  the  GAMMA 
NMR  libraries  (38)  from  a  3T  localized  2D  correlated  spectroscopy 
sequence  (39)  with  the  following  parameters:  100fq  increments, 
1024  points  in  f2,  TR/TE=  1.5  s/30  ms,  and  spectral  bandwidths  of 
1250  and  2000  Hz  along  Fq  and  F2,  respectively.  Each  2D  spectrum 
was  line  broadened  by  1 0  Hz  and  apodized  by  a  sine-squared  filter 
along  fq  and  a  skewed  sine-squared  filter  with  skew  parameter  0.5 
along  t2.  No  baseline  corrections  were  performed  on  the  spectra. 


Because  the  additive  noise  was  random,  each  SNR  was 
simulated  and  reconstructed  20  times  per  sample  mask  to 
account  for  random  fluctuations  in  the  reconstruction.  The 
sampling  masks  were  created  using  the  2D  Poisson-gap  method 
described  earlier. 

MRSI 

The  breasts  of  three  healthy  volunteers  were  scanned  using  the 
4D  EP-COSI  sequence  on  a  Siemens  (Siemens  AG,  Erlangen,  Germany) 
3T  Trio  scanner  with  the  following  parameters:  voxel  size, 
1x1x1  cm3,  50  q  increments,  TR/TE/averages  =  1 .5  s/30  ms/1 ,  field 
of  view,  16x  16cm2  FOV,  and  spectral  bandwidths  of  1250Hz  and 
1 1 90  Hz  along  Fq  and  F2,  respectively.  Each  breast  was  scanned  twice: 
a  4x  prospective  NUS  scan  and  a  fully  sampled  scan  using  the  same 
field  of  view  and  shim.  The  NUS  scan  took  5  min  to  complete  and  the 
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Figure  3.  Simulated  quad  phantom  illustration.  Top:  spatial  distribution  of  tCho  [total  choline  (choline  +  glycerophosphocholine  +  phosphocholine)], 
Cr  (creatine),  Glx  (glutamate  and  glutamine)  and  Asp  (aspartate)  diagonal  peaks  when  fully  sampled.  Bottom:  spatial  distribution  of  the  same  tCho, 
Cr,  Glx  and  Asp  diagonal  peaks  of  the  4x  non-uniform  under-sampled  (NUS)  zero-augmented  dataset. 
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fully  sampled  scan  took  20  min.  Both  scans  were  first  apodized 
using  a  sine-squared  filter  along  fq  and  a  skewed  sine-squared  filter 
with  skew  parameter  0.5  along  f2.  No  baseline  corrections  were 
performed  on  the  in  vivo  breast  data.  The  fully  sampled  scans  were 
then  retrospectively  under-sampled  4x  using  the  same  mask  as 
employed  in  the  prospective  scan  shown  in  Fig.  4,  and  both  NUS 
datasets  were  then  reconstructed  using  MaxEnt  and  CS. 


accuracy  with  respect  to  a  fully  sampled  reference  dataset  that 
increases  as  the  two  datasets  become  more  dissimilar.  The  RMSE 
was  calculated  in  the  spatial,  spectral  domain  as: 

RMSE  =  1  ^Kldatal-lfulll)2  [9] 


RESULTS 

MR  simulations 

Quantitative  results  for  the  MaxEnt  and  CS  reconstructions  of  the 
simulated  4D  EP-COSI  dataset  at  different  NUS  rates  and  SNRs 
are  shown  in  Fig.  5.  The  top  panel  shows  the  mean  RMSE  versus 
SNR  for  zero-augmented  and  reconstructed  datasets  at  each 
NUS  rate.  The  RMSE  provides  an  estimate  of  the  reconstruction 


Figure  4.  Non-uniform  under-sampling  (NUS)  mask  used  to  under-sample 
the  ky-ti  plane  4x  in  Figs  8  and  9. 
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Figure  5.  Metrics  comparing  the  zero-augmented,  maximum  entropy 
(MaxEnt)-reconstructed  and  compressed  sensing  (CS)-reconstructed 
four-dimensional  (4D)  echo-planar  correlated  spectroscopic  imaging 
(EP-COSI)  simulated  data.  Top:  overall  root-mean-square  errors  (RMSEs) 
of  each  dataset  versus  signal-to-noise  ratio  (SNR)  for  5x,  2.5x,  1.67x  and 
1.25x  non-uniform  under-sampling  (NUS)  rates.  Bottom:  tCho  [total 
choline  (choline  +  glycerophosphocholine  +  phosphocholine)],  Cr  (crea¬ 
tine),  Glx  (glutamate  and  glutamine)  and  Asp  (aspartate)  metabolite- 
specific  RMSEs  for  each  dataset  versus  SNR  for  5x,  2.5x,  1.67x  and 
1.25x  NUS  rates. 


where  N  is  the  number  of  data  points,  'full'  is  the  fully  sampled 
dataset  and  'data'  is  the  zero-augmented  or  reconstructed 
dataset.  Error  bars  are  not  shown  because  the  standard  devia¬ 
tions  were  three  to  four  orders  of  magnitude  smaller  than  the 
mean  RMSEs  and  did  not  vary  noticeably  over  the  NUS  rate  or 
SNR.  As  can  be  seen,  the  RMSE  of  the  zero-augmented  dataset 
increases  as  the  NUS  rate  increases,  but  does  not  vary  consider¬ 
ably  with  SNR,  except  at  low  NUS  rates.  Both  the  CS  and  MaxEnt 
reconstructions  show  large  decreases  in  RMSE  at  each  SNR  and 
NUS  rate,  but,  at  low  SNRs,  RMSE  begins  to  rise.  At  low  NUS  rates, 
CS  and  MaxEnt  have  comparable  RMSEs  at  each  SNR;  however, 
at  higher  NUS  rates,  the  RMSEs  for  the  MaxEnt  reconstruction 
are  lower  than  for  CS,  and  this  difference  increases  with  SNR. 

The  bottom  panel  of  Fig.  5  shows  the  average  RMSE  of  the 
diagonal  and  cross-peaks  of  tCho,  Glx,  Cr  and  Asp  versus  SNR  of 
the  zero-augmented  and  reconstructed  datasets.  The  RMSEs 
were  calculated  over  the  metabolite  peaks  at  the  ppm  locations 
listed  in  Table  1.  Each  RMSE  was  calculated  only  over  the  four 
spatially  distributed  voxels  for  each  metabolite.  For  example, 
the  tCho  RMSEs  were  calculated  over  voxels  (2,  2),  (2,  3),  (3,  2) 
and  (3,  3),  as  illustrated  in  Fig.  3.  Therefore,  these  RMSEs  reflect 
local  changes  to  the  metabolite  peak  lineshape  and  amplitude 
caused  by  the  spatial,  spectral  aliasing  along  the  ky-ti  plane  from 
the  NUS  and  reconstruction. 

All  of  the  metabolite  RMSEs  in  the  bottom  panel  of  Fig.  5  show 
similar  trends  over  SNR  as  the  overall  RMSEs  in  the  top  panel  of 
Fig.  5  at  each  NUS  rate.  The  MaxEnt  and  CS  reconstructions  have 
lower  RMSEs  than  the  zero-augmented  datasets,  indicating  that 
the  metabolite  peak  lineshapes  and  amplitudes  are  being 
properly  reconstructed.  The  MaxEnt  reconstructions  have  lower 
RMSEs  for  many  metabolites  than  the  CS  reconstructions  at 
higher  NUS  rates,  but,  at  lower  rates,  their  RMSEs  are  roughly 
equivalent. 

Figure  6  shows  a  1 D  cross-section  of  the  fully  sampled,  zero- 
augmented  and  reconstructed  spectra  for  high  and  low  SNR 
simulated  spectra  at  1.25x  and  5x  NUS,  respectively.  The  ID 
cross-section  is  indicated  by  the  broken  line  across  F2  at  F-,  =  3.65 
ppm  in  Fig.  7,  and  any  NUS  artifacts  are  from  aliased  peaks  above 
and  below  the  line,  not  peaks  shown  in  the  cross-section.  The 
high  SNR,  1.25x  NUS  zero-augmented  spectrum  shows  only 
small  deviations  from  the  fully  sampled  spectrum,  but  they  are 
clearly  visible  in  the  inset.  Both  reconstructions  restored  the 
baseline  to  the  level  of  the  fully  sampled  spectrum  and 
preserved  the  amplitude  and  lineshapes  of  the  peaks  in  the  full 
cross-sections.  The  artifacts  in  the  low  SNR,  5x  NUS  zero-augmented 
spectrum  show  significantly  reduced  peak  amplitudes,  broader 
linewidths  and  Gibbs  ringing  along  the  baseline.  Both  recon¬ 
structions  successfully  restored  the  linewidths  of  the  peaks  in 
the  full  cross-sections  and  removed  the  Gibbs  ringing  shown 
in  the  insets;  however,  MaxEnt  was  generally  better  at 
restoring  the  peak  amplitude  as  indicated  by  the  red  arrows. 
Many  of  the  real  and  imaginary  peak  amplitudes  in  the  CS 
reconstructions  were  lower  than  those  of  the  MaxEnt  recon¬ 
structions,  and  the  baseline  for  CS  was  also  slightly  lower. 
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Table  1.  Two-dimensional  peak  locations  (ppm)  for  selected  metabolites 

tCho  Glx  Cr  Asp 

Diagonals  (ppm)  (3.2, 3.2),  (3.5, 3.5),  (4.0,4.0),  (4.3,43)  (23,2.3),  (3.7, 3.7)  (3.0, 3.0),  (3.9, 3.9)  (2.8, 2.8),  (3.9, 3.9) 

Cross-peaks  (ppm)  (3.5, 4.0),  (3.5,4 .3),  (4.0,3.5),  (43,3.5)  (23,3.7),  (3.7,23)  N/A  (2.8, 3.9),  (3.9, 2.8) 

Asp,  aspartate;  Cr,  creatine;  Glx,  glutamate  and  glutamine;  N/A,  not  applicable;  tCho,  total  choline  (choline + 
glycerophosphocholine  +  phosphocholine). 


Figure  6.  One-dimensional  (ID)  magnitude,  real  and  imaginary  cross-sections  of  fully  sampled,  zero-augmented  and  reconstructed  Glx  (glutamate  and 
glutamine)  spectra.  Top:  cross-sections  from  high-signal-to-noise  ratio  (SNR)  spectra  1.25x  non-uniform  under-sampled  (NUS)  zero-augmented  and 
reconstructed  using  maximum  entropy  (MaxEnt)  or  compressed  sensing  (CS).  Bottom:  cross-sections  from  low-SNR  spectra  5x  NUS  zero-augmented 
and  reconstructed  using  MaxEnt  or  CS.  Insets:  magnified  cross-sections  of  spectral  baselines. 


Glx 


Figure  7.  Fully  sampled  two-dimensional  (2D)  correlated  spectroscopy 
(COSY)  Glx  (glutamate  and  glutamine)  spectrum,  with  the  broken  line 
across  F2  at  Fq  =  3.65  ppm  indicating  the  one-dimensional  (1 D)  cross-section 
shown  in  Fig.  6. 

NUS  of  4D  EP-COSI  in  human  breast 

The  results  from  a  prospective  4x  NUS  EP-COSI  scan  ofa  31-year-old 
healthy  human  breast  that  was  MaxEnt-  and  CS-reconstructed  are 


shown  in  Fig.  8A,  B,  respectively,  with  the  zero-augmented  data 
shown  in  Fig.  8C.  The  contour  levels  employed  in  the  fully  sampled 
results  in  Fig.  9A  were  used  in  Fig.  8.  The  mask  used  to  under-sample 
the  /cy-fq  plane  is  shown  in  Fig.  4,  together  with  the  signal  envelopes 
for  each  dimension,  and  was  generated  using  the  2D  Poisson-gap 
method  described  earlier  in  this  article. 

Figure  8A1,  B1  shows  2D  correlated  spectroscopy  (COSY) 
spectra  extracted  from  the  MaxEnt  and  CS  reconstructions, 
respectively.  They  were  taken  from  the  fatty  breast  regions 
highlighted  in  Fig.  8A2,  B2.  They  clearly  show  the  lipid  diagonal 
peaks,  olefinic  fat  (UFD),  methyl  fat  (FMETD)  and  fat  (FAT/FAT2/ 
FAT3),  and  the  cross-peaks,  unsaturated  fatty  acid  right  (UFR), 
unsaturated  fatty  acid  left  (UFL),  and  triglyceryl  fat  (TGFR)  (5).  The 
spatial  distribution  of  the  UFL/UFR  cross-peaks  from  the  recon¬ 
structions  is  shown  in  Fig.  8A2,  B2  with  the  MaxEnt  reconstruction's 
spatial  distribution  overlaid  on  the  anatomical  MR  image. 

Figure  8C1  shows  the  same  2D  COSY  spectrum  as  in  Fig.  8A1,  B1 
with  4x  NUS  applied  to  the  ky-t q  plane  using  the  mask  in  Fig.  4; 
however,  no  MaxEnt  or  CS  reconstruction  was  used.  The  spatial, 
spectral  incoherent  artifacts  from  NUS  manifest  as  smeared  peaks 
along  Fq,  which  is  illustrated  by  the  collapse  of  the  peaks  in  the 
ID  projection  of  the  Fq  dimension  on  the  right.  The  aliasing  of  the 
large  diagonal  fat  peaks  around  (F2  =  2ppm,  Fq=2ppm)  obscures 
the  much  smaller  UFL/UFR  cross-peaks  around  (F2  =  2.1  ppm, 
Fq  =  5.4  ppm).  Figure  8C2  shows  the  spatial  distribution  of  the 
UFL/UFR  cross-peaks  and  how  spatial  artifacts  from  the  under¬ 
sampling  of  ky-fq  manifest  as  errant  peaks  in  adjacent  voxels. 
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Figure  8.  Prospective  four-dimensional  (4D)  echo-planar  correlated  spectroscopic  imaging  (EP-COSI)  results.  Top:  selected  two-dimensional  (2D) 
correlated  spectroscopy  (COSY)  spectrum  from  a  4D  EP-COSI  scan  of  healthy,  fatty  breast  highlighted  in  the  bottom  images  for  maximum  entropy 
(MaxEnt)-reconstructed  (A1),  compressed  sensing  (CS)-reconstructed  (B1),  and  4x  non-uniform  under-sampled  (NUS)  zero-augmented  (Cl)  data. 
Bottom:  spatial  distribution  of  the  unsaturated  fatty  acid  left/unsaturated  fatty  acid  right  (UFL/UFR)  cross-peaks  highlighted  in  the  2D  COSY  spectrum 
for  MaxEnt-reconstructed  (A2),  CS-reconstructed  (B2),  and  4x  NUS  zero-augmented  (C2)  data.  FAT/FAT2/FAT3,  fat;  FMETD,  methyl  fat;  TGFR,  triglyceryl 
fat;  UFD,  olefinic  fat. 
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Figure  9.  Retrospective  four-dimensional  (4D)  echo-planar  correlated  spectroscopic  imaging  (EP-COSI)  results.  Top:  selected  two-dimensional  (2D) 
correlated  spectroscopy  (COSY)  spectrum  from  a  4D  EP-COSI  scan  of  healthy,  fatty  breast  highlighted  in  the  bottom  images  for  fully  sampled  (A1), 
maximum  entropy  (MaxEnt)-reconstructed  (B1),  and  compressed  sensing  (CS)-reconstructed  (Cl)  data.  Bottom:  spatial  distribution  of  the  unsaturated 
fatty  acid  left/unsaturated  fatty  acid  right  (UFL/UFR)  cross-peaks  highlighted  in  the  2D  COSY  spectrum  for  fully  sampled  (A2),  MaxEnt-reconstructed 
(B2),  and  CS-reconstructed  (C2)  data.  FAT/FAT2/FAT3,  fat;  FMETD,  methyl  fat;  TGFR,  triglyceryl  fat;  UFD,  olefinic  fat. 


Comparing  the  MaxEnt-  and  CS-reconstructed  spectra  in 
Fig.  8A1,  B1,  all  of  the  significant  diagonals  and  cross-peaks  are 
fully  resolved  in  both  datasets  with  qualitatively  similar 
linewidths  and  amplitudes.  The  only  major  differences  between 
them  are  that  the  amplitudes  of  the  f,  ridges  centered  at  (F2  =  1  -3 
ppm,  F-,  =  1.3ppm)  and  (F2=1.0ppm,  F^I.Oppm)  for  CS  are 
lower  than  for  MaxEnt. 

The  results  from  a  fully  sampled  EP-COSI  scan  of  the  same 
healthy  breast  as  shown  in  Fig.  8,  which  was  retrospectively  4x 
under-sampled  using  the  same  mask  and  MaxEnt  and  CS 
reconstructions,  are  illustrated  in  Fig.  9.  Figure  9A  shows  the  fully 
sampled  data  and  Fig.  9B,  C  shows  the  MaxEnt  and  CS  recon¬ 
structions,  respectively.  As  the  same  field  of  view  was  used  for 
both  scans,  the  spectra  at  the  top  of  Fig.  9  show  the  same  2D 
COSY  spectra  as  those  at  the  top  of  Fig.  8,  and  the  bottom 
shows  the  spatial  distribution  of  the  UFL/UFR  cross-peaks.  The 
same  contour  levels  used  for  the  fully  sampled  data  in  Fig.  9A 
were  also  employed  for  the  MaxEnt-  and  CS-reconstructed 
results  in  Fig.  9B,  C. 

NMR  Biomed.  2014;  27:  191-201 


As  can  be  seen,  all  of  the  peaks  in  the  fully  sampled  spectrum 
in  Fig.  9A1  are  completely  resolved  in  both  the  prospective  and 
retrospective  reconstructions  shown  in  Figs  8B1,  Cl,  9B1,  Cl; 
their  positions,  linewidths,  amplitudes  and  spectral  resolutions 
are  all  qualitatively  comparable.  However,  the  spatial  distribu¬ 
tions  of  the  fully  sampled  and  retrospective  NUS  results  show 
better  agreement  than  the  prospective  NUS  results  in  the  upper 
region  of  the  breast;  the  excited  volume  of  the  prospective  NUS 
results  is  one  row  smaller  than  that  of  the  retrospective  NUS 
results,  as  indicated  by  the  arrows  in  Figs  8  and  9.  This  change 
was  also  observed  in  the  non-water-suppressed  scans  taken 
prior  to  the  prospective  NUS  and  fully  sampled  scans,  and  there¬ 
fore  cannot  be  an  artifact  of  the  reconstruction. 

Table  2  shows  the  F-,  full  width  at  half-maximum  (FWHM)  and 
peak  amplitudes  of  the  zero-augmented  and  reconstructed 
dataset  magnitude  peaks  in  Figs  8  and  9.  For  comparison,  they 
are  normalized  by  the  fully  sampled  peak  amplitudes  and  FWHM 
from  Fig.  9,  so  that  values  greater  than  unity  are  larger  than  the 
fully  sampled  value.  As  can  be  seen,  there  are  quantitative 
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Table  2.  Relative  full  width  at  half-maximum  (FWFIM)  along  Fq  and  amplitude  of  metabolite  peaks  for  zero-augmented,  maximum 
entropy  (MaxEnt)-reconstructed  and  compressed  sensing  (CS)-reconstructed  data  from  the  voxel  shown  in  Figs  8  and  9.  Values  are 
normalized  by  the  fully  sampled  peak  FWFIM  and  amplitudes 

Prospective  amplitude 

FAT 

FAT2 

FAT3 

FMETD 

TGFR 

UFD 

UFL 

UFR 

MaxEnt 

1.249 

1.342 

1.232 

1.614 

0.951 

1.125 

1.060 

1.078 

CS 

1.427 

1.372 

1.263 

1.665 

0.982 

1.142 

1.150 

1.161 

Prospective  FWFIM 

MaxEnt 

0.818 

0.909 

0.818 

1.000 

0.818 

1.000 

0.900 

0.727 

CS 

0.727 

0.909 

0.818 

0.900 

0.909 

1.000 

0.900 

0.727 

Retrospective  amplitude 

Zero-augmented 

0.634 

0.701 

0.662 

0.639 

0.639 

0.627 

0.518 

0.845 

MaxEnt 

1.061 

1.026 

1.090 

1.010 

1.019 

1.029 

0.942 

0.989 

CS 

1.237 

1.138 

1.017 

0.961 

1.077 

0.995 

1.071 

1.043 

Retrospective  FWFIM 

Zero-augmented 

1.273 

1.273 

1.182 

1.30 

1.00 

1.30 

1.40 

1.010 

MaxEnt 

0.818 

0.909 

0.818 

1.00 

0.818 

1.00 

1.00 

0.818 

CS 

0.727 

0.818 

1.010 

1.00 

0.727 

1.10 

0.80 

0.727 

FAT/FAT2/FAT3,  fat;  FMETD,  methyl  fat;  TGFR,  triglyceryl  fat;  UFD,  olefinic  fat;  UFL,  unsaturated  fatty  acid  left;  UFR,  unsaturated  fatty  acid  right. 

differences  between  the  CS-  and  MaxEnt-reconstructed  peak 
lineshapes.  The  zero-augmented  dataset  has  broader,  shorter 
peaks  as  expected,  and  both  the  CS-  and  MaxEnt-reconstructed 
peak  lineshapes  are  closer  to  the  fully  sampled  data;  however, 
the  MaxEnt  peak  lineshapes  are  almost  all  closer  to  the  fully 
sampled  data  than  are  the  CS  peak  lineshapes,  which  are 
narrower  and  taller  than  both  the  MaxEnt  and  fully  sampled 
data.  The  increase  in  peak  amplitudes  in  the  CS  reconstruction 
is  greater  in  the  larger  peaks  (FAT,  FAT2  and  FAT3)  than  in  the 
smaller  peaks  (UFD,  UFL,  UFR)  when  compared  with  MaxEnt. 

Additional  quantitative  differences  between  the  MaxEnt  and 
CS  reconstructions  are  illustrated  in  Table  3,  which  shows  the 
range  of  values  for  the  (UFL+ UFR)/(FAT3  +  FAT2)  integrated 
peak  area  ratios  and  average  errors  from  the  fully  sampled  ratios 
for  the  central  6x6  voxels  of  the  three  healthy  human  breasts. 
As  can  be  seen,  NUS  caused  the  ratios  to  vary  considerably  from 
the  fully  sampled  ratios  with  a  high  mean  ratio  error.  The 
reconstructed  ratios  show  much  better  agreement  with  the  fully 


sampled  ratios  and  have  a  much  smaller  mean  ratio  error  in  all 
three  scans  than  do  the  zero-augmented  data.  Flowever,  the 
mean  ratio  errors  for  the  MaxEnt  reconstruction  are  almost  all 
smaller  than  their  CS  counterparts,  indicating  that  the  peak  ratio 
was  more  accurately  reconstructed  by  MaxEnt.  The  retrospective 
reconstruction  results  are  only  slightly  better  than  the  prospec¬ 
tive  results,  which  is  not  surprising,  given  the  excitation  volume 
differences  seen  in  Figs  8  and  9. 


DISCUSSION 

The  simulated  4D  EP-COSI  dataset  results  in  Fig.  5  demonstrate 
how  the  MaxEnt  and  CS  reconstructions  perform  at  different 
SNRs  and  NUS  rates.  Both  reconstructions  decrease  the  RMSE 
significantly  at  each  SNR  and  NUS  rate,  compared  with  the 
zero-augmented  dataset,  but  begin  to  increase  at  very  low 
SNR,  which  indicates  that  they  are  unable  to  fully  reconstruct 


Table  3.  Range  and  mean  error  of  (UFL  + UFR)/(FAT3  +  FAT2)  integrated  peak  area  ratios  for  fully  sampled,  zero-augmented, 
maximum  entropy  (MaxEnt)-reconstructed,  and  compressed  sensing  (CS)-reconstructed  data  from  three  healthy  breasts 


Ratio  range 

Mean  ratio  error 

Ratio  range 

Mean  ratio  error 

Ratio  range 

Mean  ratio  error 

Breast  1 

Breast  2 

Breast  3 

Full 

Zero-augmented 
Retrospective  MaxEnt 
Prospective  MaxEnt 
Retrospective  CS 
Prospective  CS 

0.1277-0.0456 

0.2162-0.1496 

0.1061-0.0281 

0.1049-0.0216 

0.1103-0.0275 

0.1084-0.0211 

N/A 

0.0823  ±0.01 98 

0.0244  ±0.0071 

0.0326  ±0.0083 

0.0246  ±  0.0073 

0.0334  ±0.0082 

0.1984-0.0656 

0.4699-0.3082 

0.1576-0.0413 

0.1215-0.0420 

0.1781-0.1004 

0.1276-0.0425 

N/A 

0.2640  ±  0.0734 

0.0254  ±0.01 42 

0.041 5  ±0.01 64 

0.0357  ±0.0235 

0.0397  ±0.01 40 

0.3542-0.0277 

1.0077-1.1785 

0.4134-0.0291 

0.3422-0.0414 

0.4114-0.0269 

0.3594-0.0374 

N/A 

0.2047  ±0.091 6 

0.021 5  ±0.01 49 

0.0404  ±0.0463 

0.0268  ±0.01 75 

0.0434  ±0.0487 

FAT2/FAT3,  fat;  N/A,  not  applicable;  UFL,  unsaturated  fatty  acid  left;  UFR,  unsaturated  fatty  acid  right. 
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the  data  when  features  and  sampling  artifacts  are  obscured  by 
high  levels  of  noise.  The  data  fidelity  constraint  in  Equations  [1] 
and  [2]  determines  how  closely  the  reconstructed  points  must 
be  to  the  sampled  points  within  the  standard  deviation  of  noise, 
which  increases  as  the  noise  floor  increases.  This  increase  in  the 
noise  floor  effectively  'loosens'  the  fidelity  constraint,  which  al¬ 
lows  the  reconstructed  points  to  deviate  from  their  sampled  coun¬ 
terparts  and  increases  the  entropy  or  sparsity  of  the  reconstructed 
spectrum  by  narrowing  the  peak  linewidths  and  de-noising  smaller 
features.  This,  in  turn,  increases  the  non-linearity  and  RMSE  of  the 
reconstruction  because  of  the  loose  fidelity  constraint  (41).  The 
data  fidelity  constraint  can  be  'tightened'  beyond  the  standard  de¬ 
viation  of  the  noise  in  an  effort  to  reduce  the  RMSE  and  reconstruc¬ 
tion  non-linearity,  but  this  prevents  the  Cambridge  algorithm  from 
completely  removing  the  spatial,  spectral  NUS  artifacts  close  to  the 
noise  floor,  which  could  potentially  obscure  small  features  (21). 

The  metabolite  RMSEs  in  Fig.  5  show  that  the  reconstructions 
offer  significant  improvements  in  the  amplitude  and  lineshape  of 
individual  peaks  over  the  zero-augmented  data,  even  without 
using  methods  to  reduce  the  reconstruction  non-linearity  (32). 
CS  and  MaxEnt  produce  very  similar  results  at  low  NUS  rates, 
but  MaxEnt  generally  has  a  lower  RMSE  at  higher  NUS  rates  than 
this  implementation  of  CS.  As  demonstrated  by  Figs  6  and  7, 
both  reconstruction  methods  successfully  restored  the  peaks 
and  removed  the  Gibbs  ringing  artifacts;  however,  CS  had  a 
tendency  to  over-smooth  smaller  features  close  to  the  noise 
floor  and  narrow  peak  linewidths.  It  was  this  over-smoothing 
that  reduced  the  amplitude  of  the  CS-reconstructed  fat  tails  at 
(F2=1.3ppm,  Fq  =  1.3ppm)  and  (F2=1.0ppm,  Fq  =  1.0ppm)  in 
Fig.  8.  The  differences  in  peak  linewidths  between  MaxEnt  and 
CS  contributed  to  the  disparity  in  RMSE  values  at  higher  NUS 
rates,  but  the  main  contributing  factor  was  the  change  in  the 
noise  floor  in  the  non-peak  regions  of  the  spectra.  The  noise  floor 
was  slightly  reduced  over  the  entire  4D  dataset  in  the  CS  recon¬ 
structions,  which  was  the  vast  majority  of  points  in  the  volume; 
therefore,  small  changes  in  the  noise  floor  had  large  effects  on 
RMSE.  The  over-smoothing  in  the  CS  reconstructions  may  have 
been  caused  by  choosing  e  =  def! 1000  in  Equation  [4],  and  will 
be  investigated  in  the  future,  although  previous  work  has  used 
similar  values  with  success  (13). 

The  healthy  human  breast  results  in  Figs  8  and  9  show  that  CS 
and  MaxEnt  reconstructions  work  well  for  prospective  and 
retrospective  NUS,  filtered,  in  vivo  EP-COSI  scans.  There  were 
minor  differences  between  the  prospective  and  retrospective 
reconstructions;  however,  these  can  be  attributed  to  intra-scan 
variations  in  the  excitation  volumes,  as  these  differences  were 
reflected  in  the  MaxEnt  and  CS  reconstructions,  as  well  as  in 
the  non-water-suppressed  scans.  The  greater  SNR  loss  caused 
by  T2*  decay  which  was  not  present  in  the  simulated  4D  EP-COSI 
dataset  did  not  reduce  the  efficacy  of  the  reconstructions.  There 
was  still  sufficient  SNR  in  the  time  domain  to  reconstruct  the 
in  vivo  diagonal  and  cross-peaks  in  the  spectral  domain. 

Although  the  MaxEnt  and  CS  reconstructions  of  healthy  human 
breast  were  qualitatively  similar,  Tables  2  and  3  illustrate  quanti¬ 
tative  differences  between  them.  Both  reconstruction  methods 
improved  the  ratios  of  the  lipid  peaks  in  Table  3  and  the  ampli¬ 
tudes  and  FWHM  of  major  peaks  in  Table  2,  which  shows  that 
they  were  able  to  reconstruct  the  large,  aliased  diagonals,  as  well 
as  the  smaller  cross-peaks  that  were  obscured  by  the  diagonals 
aliasing  over  the  ky-fq  plane.  However,  as  shown  in  Table  2,  the 
CS-reconstructed  peaks  were  narrower  with  higher  amplitudes 
than  their  MaxEnt  counterparts.  This  discrepancy  increased  with 


peak  amplitude,  indicating  a  higher  degree  of  non-linearity  in 
the  CS  reconstruction  relative  to  MaxEnt.  Using  the  relative  peak 
amplitudes  and  FWHM  values  from  Table  2  to  calculate  the  peak 
area,  instead  of  the  integrated  peakarea  that  was  used  forTable  3, 
the  relative  (UFL+ UFR)/(FAT3  + FAT2)  ratios  for  MaxEnt  and  CS 
for  the  retrospective  reconstructions  are  0.960  and  0.825,  respec¬ 
tively.  Because  the  FAT3  and  FAT2  peaks  from  the  CS  reconstruc¬ 
tion  are  relatively  larger  than  the  MaxEnt  peaks  from  the 
increased  non-linearity,  their  relative  ratio  with  the  UFL  and  UFR 
peaks  is  smaller.  This  increased  non-linearity  was  a  contributing 
factor  in  the  larger  mean  ratio  errors  in  Table  3  for  the  CS  recon¬ 
structions  relative  to  the  MaxEnt  reconstructions. 

The  under-sampled  data  in  Fig.  8C1,  C2  show  artifacts  spread 
along  Fq  and  Y,  as  well  as  reduced  spectral  resolution  along  Fq 
and  larger  FWHMs  in  Table  2,  caused  by  convolution  with  the 
broad  NUS  PSF.  The  homogeneous  nature  of  healthy  fatty  breast 
spectra,  coupled  with  the  inherently  low  spatial  resolution  of  4D 
EP-COSI,  made  it  difficult  to  determine  from  the  figures  whether 
the  spatial  resolution  along  Y  decreased  as  a  result  of  the  NUS 
PSF.  However,  it  is  clear  that  the  effects  of  the  NUS  PSF  along 
Fq  were  removed  by  MaxEnt  and  CS  by  the  narrower  FWHM 
values  in  Table  2  and  the  lack  of  NUS  artifacts  in  Figs  8A1,  Bl, 
9A1 ,  B1 .  The  errant  spectral  peaks  in  the  spatial  distribution  were 
removed  in  Figs  8A2,  B2,  9A2,  B2,  suggesting  that  the  spatial  PSF 
along  Y  was  also  improved.  Any  spectral  bleed  from  the  spatial 
PSF  of  the  EP-COSI  pulse  sequence  along  X  was  orthogonal  to 
the  effect  of  the  NUS  PSF  along  Y,  and  was  not  affected  by  the 
MaxEnt  and  CS  reconstructions. 

The  results  in  Figs  5-9  indicate  that  MaxEnt  and  CS  produced 
qualitatively  similar  reconstructions;  however,  the  MaxEnt  results 
were  quantitatively  better  by  a  small  margin  as  discussed  above. 
This  is  not  surprising,  given  that  the  objective  functions  of 
CS  and  MaxEnt  are  similar,  but  there  are  minor  differences 
between  them  as  shown  in  Fig.  1.  CS  uses  the  lrnorm  of  the 
reconstructed  spectrum  in  some  transform  domain,  and  MaxEnt 
uses  entropy,  which  is  a  log-sum  function  that  can  be  rewritten 
as  a  reweighted  lrnorm,  I  wr  |m,-|,  where  w,-  is  log[f(m,)],  the  log 
of  a  function  of  the  reconstructed  spectrum.  Previous  work  has 
shown  that  reweighted  lrnorm  objective  functions  can  outperform 
/j-norm-based  CS  reconstruction  (42),  and  direct  comparisons  bet¬ 
ween  MaxEnt  and  lrnorm- based  CS  reconstruction  have  shown 
them  to  be  qualitatively  equivalent  (43).  This  is  the  first  known  work 
to  show  quantitative  comparisons  between  these  techniques,  how¬ 
ever,  and  further  research  into  their  relative  performance  is  ongoing. 

The  Cambridge  algorithm  used  to  solve  the  MaxEnt  and  CS 
reconstruction  problems  was  demonstrated  to  be  robust  against 
different  levels  of  SNR  and  NUS  rates  for  the  simulated  and 
in  vivo  datasets.  There  are  other  irnorm  solvers  available  for 
the  CS  reconstruction  problem;  however,  many  require  parame¬ 
ter  tuning  for  different  datasets  in  order  to  find  the  optimal 
reconstruction  parameters  (44,45).  The  Cambridge  algorithm 
does  not  have  any  tuning  parameters  that  must  be  adjusted  to 
find  the  optimal  reconstruction  for  a  dataset,  which  offers  a  sub¬ 
stantial  workflow  improvement  over  other  solvers.  Although  the 
Cambridge  algorithm  can  be  modified  to  solve  the  CS  lrnorm 
reconstruction  problem,  it  takes,  on  average,  5-10  times  longer 
to  converge  as  the  MaxEnt  problem,  which  took  7-10  min,  on 
average,  using  a  64-bit  dual-core,  3.4-GHz  Core  i7  processor  with 
16GB  RAM.  Therefore,  we  do  not  recommend  its  use  as  an 
lrnorm  solver;  it  was  only  used  for  the  current  work  in  order  to 
compare  results  for  MaxEnt  and  CS  reconstruction  that  were 
not  biased  by  different  solver  implementations.  However, 
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because  it  is  relatively  fast,  robust  to  SNR  changes,  does  not 
require  parameter  tuning,  and  provides  MaxEnt  results  that  were 
equal  to  or  better  than  those  of  CS,  the  Cambridge  algorithm  is 
well  suited  as  a  MaxEnt  solver. 

The  Poisson-gap  sampling  masks  used  in  these  experiments 
were  generated  by  a  random  Poisson  distribution,  which 
injects  a  degree  of  uncertainty  into  the  reconstruction.  It  has 
proven  to  be  a  reliable  technique  that  generates  masks  with 
desirable  PSFs,  as  shown  in  Fig.  2,  and  the  RMSEs  of  the 
reconstructed  datasets  using  different  Poisson-gap  sample 
masks  are  stable  (28).  Previous  attempts  by  our  group  to  use 
deterministic  masks  that  were  not  randomly  generated  were 
difficult  to  optimize  and  suffered  from  coherent  aliasing, 
which  cannot  be  removed  by  MaxEnt  or  CS  reconstruction 
(32,46).  Recent  work  in  NUS  multi-dimensional  NMR  datasets 
using  deterministic  sample  masks  has  shown  promise  and 
could  be  adapted  to  4D  MRSI  (47). 

Because  of  the  random  nature  of  the  Poisson-gap  sampling 
masks,  they  were  chosen  by  an  empirical  heuristic  that  minimized 
the  width  of  the  central  peak,  the  total  power  of  the  incoherent 
artifacts,  and  the  ratio  of  the  largest  artifact  peak  to  the  central  peak 
in  the  mask  PSF.  They  followed  a  skewed-sine  bell  modulation  func¬ 
tion  to  maximize  the  reconstruction  SNR,  but  were  not  optimized 
for  specific  metabolites  or  post-processing  spectral  filters.  In  our 
experiments,  we  observed  that  mismatches  between  a  sample 
modulation  function  and  the  filtered  signal  envelope  of  a  metabo¬ 
lite  results  in  failure  to  sufficiently  sample  the  high-SNR  points  along 
f,  and  prevents  the  full  metabolite  peak  area  from  being 
reconstructed  (31,32,48).  Therefore,  it  should  be  emphasized  that 
the  sample  mask  is  crucial  to  the  SNR  of  each  reconstructed  metab¬ 
olite  for  4D  EP-COSI  data,  and  there  is  a  dependence  on  the  shape 
of  standard  spectral  filters  applied  prior  to  reconstruction  (33). 

Further  research  into  reducing  the  non-linearity  of  the 
reconstructed  peaks  to  make  accurate  quantification  possible, 
and  comparisons  with  additional  CS  methods,  is  ongoing.  In  addi¬ 
tion,  future  papers  will  address  the  use  of  Poisson-gap  versus 
deterministic  sample  masks,  and  the  optimization  of  the  modula¬ 
tion  functions  for  specific  metabolites  and  different  spectral  filters. 

CONCLUSIONS 

This  work  has  demonstrated  that  MaxEnt  is  a  viable  alternative  to  /,- 
norm-based  CS  reconstruction  for  accelerating  the  acquisition  of  4D 
EP-COSI  data  in  vivo.  MaxEnt  provided  reconstructions  equal  to  or  bet¬ 
ter  than  those  of  CS,  and  the  robust  nature  of  the  Cambridge  algorithm 
without  the  need  for  parameter  tuning  makes  it  a  good  candidate  for 
clinical  use.  The  CS  and  MaxEnt  reconstructions  throughout  this 
work  were  qualitatively  similar;  however,  the  quantitative  results 
indicated  increased  non-linearity  in  the  CS  reconstruction  when 
compared  to  MaxEnt.  Simulated  4D  EP-COSI  data  provided  a  quan¬ 
titative  characterization  of  both  reconstruction  methods  at  differ¬ 
ent  NUS  rates  and  SNRs,  and  the  4x  NUS  in  vivo  EP-COSI  breast 
data  showed  that  a  clinically  viable  5-min  breast  scan  is  possible. 

Acknowledgements 

The  authors  acknowledge  the  scientific  support  of  Dr  Rajakumar 
Nagarajan  for  recording  the  in  vivo  NUS  EP-COSI  data.  This  work 
was  supported  by  a  National  Institutes  of  FHealth  (NIFH)  training 
grant  #5T15  LM07356  (BB)  and  an  IDEA  Expansion  grant  from 
the  US  Army  Department  of  Defense  (DOD)  Breast  Cancer 
Research  Program  (BCRP)#W81XWH-1 0-1 -0743  (MAT). 


REFERENCES 

1 .  Thomas  MA,  Lipnick  S,  Velan  SS,  Liu  X,  Banakar  S,  Binesh  N,  Ramadan  S, 
Ambrosio  A,  Raylman  RR,  Sayre  J,  DeBruhl  N,  Bassett  L.  Investigation  of 
breast  cancer  using  two-dimensional  MRS.  NMR  Biomed.  2008;  22(1): 
77-91. 

2.  Begley  JK,  Redpath  TW,  Bolan  PJ,  Gilbert  FJ.  In  vivo  proton  magnetic 
resonance  spectroscopy  of  breast  cancer:  a  review  of  the  literature. 
Breast  Cancer  Res.  2012;  14(2):  207. 

3.  Glunde  E,  Bhujwalla  ZM,  Ronen  SM.  Choline  metabolism  in  malig¬ 
nant  transformation.  Nat  Rev  Cancer,  2011;  11  (2):  835-848. 

4.  Aue  W,  Bartholdi  E,  Ernst  R. Two-dimensional  spectroscopy.  Application 
to  nuclear  magnetic  resonance.  J.  Chem.  Phys.  1976;  64:  2229-2246. 

5.  Lipnick  S,  Liu  X,  Sayre  J,  Bassett  LW,  DeBruhl  N,  Thomas  MA.  Combined 
DCE-MRI  and  single-voxel  2D  MRS  for  differentiation  between  benign 
and  malignant  breast  lesions.  NMR  Biomed.  2010;  23(8):  922-930. 

6.  Mansfield  P.  Spatial  mapping  of  the  chemical  shift  in  NMR.  J  Phys  D 
Appl  Phys  1983;  16:  235-238. 

7.  Mulkern  R,  Panych  L.  Echo  planar  spectroscopic  imaging.  Concepts 
Magn.  Reson.  2001;  13:  213-237. 

8.  Posse  S,  DeCarli  C,  Le-Bihan  D.  3D  echo  planar  spectroscopic  imaging 
at  short  echo  times  in  human  brain.  Radiology  1994;  192:  733-738. 

9.  Lipnick  S,  Verma  G,  Saadallah  R,  Jon  F,  Thomas  MA.  Echo  planar  cor¬ 
related  spectroscopic  imaging:  implementation  and  pilot  evaluation 
in  human  calf  in  vivo.  Magn  Reson  Med  2010;  64(4):  947-956. 

1 0.  Furuyama  JK,  Wilson  NW,  Burns  BL,  Rajakumar  N,  Margolis  D,  Thomas  MA. 
Application  of  compressed  sensing  to  multidimensional  spectroscopic 
imaging  in  human  prostate.  Magn  Reson  Med  2012;  67: 1499-1505. 

11.  Donoho  D.  Compressed  sensing.  IEEE  Trans.  Info.  Theor.  2006;  52(4): 
1289-1306. 

12.  Candes  EJ,  Romberg  J,  Tao  T.  Robust  uncertainty  principles:  exact 
signal  reconstruction  from  highly  incomplete  frequency  information. 
IEEE  Trans.  Info.  Theor.  2006;  52(2):  489-509. 

13.  Lustig  M,  Donoho  D,  Pauly  JM.  Sparse  MRI:  the  application  of  compressed 
sensing  for  rapid  MR  imaging.  Magn  Reson  Med  2007;  53(6):  1182-1195. 

14.  Block  KT,  Uecker  M,  Frahm  J.  Undersampled  radial  MRI  with  multiple 
coils.  Iterative  image  reconstruction  using  a  total  variation  constraint. 
Magn  Reson  Med  2007;  57(6):  1086-1098. 

15.  Flu  S,  Lustig  M,  Chen  A,  Crane  J,  Kerr  A,  Kelley  DAC,  Flurd  R, 
Kurhanewicz  J,  Nelson  SJ,  Pauly  JM,  Vigneron  DB.  Compressed  sensing 
for  resolution  enhancement  of  hyperpolarized  Cl  3  flyback  3D-MRSI. 
J  Magn  Reson  2008;  192:  258-264. 

16.  Gamper  U,  Boesiger  P,  Kozerke  S.  Compressed  sensing  in  dynamic 
MRI.  Magn  Reson  Med  2008;  59:  365-373. 

17.  Jung  H,  Kyunghyun  S,  Nayak  KS,  Kim  E,  Ye  JC.  k-t  FOCUSS:  a  general 
compressed  sensing  framework  for  high  resolution  dynamic  MRI. 
Magn  Reson  Med  2009;  61(1):  103-116. 

18.  Frieden  RB.  Restoring  with  maximum  likelihood  and  maximum 
entropy.  J  Opt  Soc  Am  1972;  62(4):  511-518. 

19.  Skilling  J,  Bryan  RK.  Maximum  entropy  image  reconstruction:  general 
algorithm.  Monthly  Notices  R.  Astronom.  Soc.  1 984;  21 1  (1 ):  1 1 1  -1 24. 

20.  Daniell  GJ,  Flore  PJ.  Maximum  entropy  and  NMR-a  new  approach.  J 
Magn  Reson.  1989;  4(3):  515-536. 

21.  Floch  J,  Stern  AS.  NMR  Data  Processing.  Wiley:  New  York,  NY;  1996. 

22.  Jaynes  ET.  Information  theory  and  statistical  mechanics.  Phys  Rev 
1957;  1 06(4):  620-630. 

23.  FHolland  D,  Bostock  M,  Gladden  L,  Nietlispach  D.  Fast 
multidimensional  NMR  spectroscopy  using  compressed  sensing. 
Angew  Chem  Int  Ed  2011;  50(29):  6548-6551. 

24.  Boyd  S,  Vandenberghe  L.  Convex  Optimization.  Cambridge  University 
Press:  Cambridge,  2004. 

25.  Shannon  CE.  A  mathematical  theory  of  communication.  Bell  Sys. 
Tech.  J.  1948;  27:  379-423,  623-656. 

26.  Frieden  R.  Unified  theory  for  estimating  frequency-of-occurrance 
laws  and  optical  objects.  J  Opt  Soc  Am  1983;  73(7):  927-938. 

27.  Constable  RT,  Henkelmen  RM.  Why  MEM  does  not  work  in  MR  image 
reconstuction.  Magn  Reson  Med  1990;  14:  12-25. 

28.  Flyberts  SG,  Takeuchi  K,  Wagner  G.  Poisson-gap  sampling  and  for¬ 
ward  maximum  entropy  reconstruction  for  enhancing  the  resolution 
and  senstivity  of  protein  NMR.  J  Am  Chem  Soc  2010;  132:  2145-2147. 

29.  Maciejewski  MW,  Mobli  M,  Schuyler  A,  Stern  A,  Floch  J.  Data 
sampling  in  multidimensional  NMR:  fundamentals  and  strategies. 
Top  Curr  Chem  2011;  316:  49-77. 

30.  Silverman  H,  Shapiro  R.  Alias  free  sampling  of  random  noise.  SIAM  J. 
1960;  8(2):  225-248. 


wileyonlinelibrary.com/journal/nbm 


Copyright  ©  2013  John  Wiley  &  Sons,  Ltd. 


NMR  Biomed.  2014;  27:  191-201 


NON-LINEAR  RECONSTRUCTION  METHODS  FOR  UNDER-SAMPLED  MRSI  DATA  IN  VIVO 


NMR 

IN  BIOMEDICINE 


31.  Mobli  M,  Hoch  JC.  Maximum  entropy  spectral  reconstruction  of 
nonuniformly  sampled  data.  Concepts  Magn.  Reson.  2008;  32A(6): 
436-448. 

32.  Schmieder  P,  Stern  A,  Wagner  G,  Hoch  J.  Applications  of  nonlinear  sam¬ 
pling  schemes  to  COSY-type  spectra.  J  Biomol  NMR  1993;  3(5):  569-576. 

33.  Ziegler  A,  Izquierdo  M,  Decorps  M.  Optimization  of  homonuclear 
two-dimensional  correlation  methods  for  in  vivo  and  ex  vivo  NMR. 
J.  Magn.  Reson.  Ser.  B  1995;  107:  10-18. 

34.  Marseille  GJ,  De  Beer  R,  Fuderer  M,  Mehlkopf  AF,  van  Ormondt  D. 
Nonuniform  phase-encode  distributions  for  MRI  scan  time  reduction. 
J.  Magn.  Reson.  Ser.  B  1996;  111:  70-75. 

35.  Devroye  L.  Non-Uniform  Random  Variate  Generation.  Springer-Verlag: 
New  York,  NY;  1 986. 

36.  Dippe  M,  Wold  E.  Antialiasing  through  stochastic  sampling.  ACM 
SIGGRAPH  Comp.  Graphics  1985;  19(3):  69-78. 

37.  Cook  R.  Stochastic  sampling  in  computer  graphics.  ACM  Trans. 
Comp.  Graphics  1986;  5(1):  51-72. 

38.  Smith  SA,  Levante  TO,  Meier  BH,  Ernst  RR.  Computer  simulations  in 
magnetic  resonance,  an  object-oriented  programming  approach.  J 
Magn  Reson  1994;  106(1):  75-105. 

39.  Thomas  MA,  Yue  K,  Binesh  N,  Davanzo  P,  Kumar  A,  Siegel  B,  Frye  M, 
Curran  J,  Lufkin  R,  Martin  P,  Guze  B.  Localized  two-dimensional  shift 
correlated  MR  spectroscopy  of  human  brain.  Magn  Reson  Med  2001; 
46(1):  58-67. 


40.  Pierpaoli  C,  Basser  PJ.  Toward  a  quantitative  assesment  of  diffusion 
anisotropy.  Magn  Reson  Med  1996;  36:  893-906. 

41.  Hoch  J,  Stern  A,  Donoho  D,  Johnstone  IM.  Maximum  entropy 
reconstruction  of  complex  (phase-sensitive)  spectra.  J  Magn  Reson 
1990;  86:  236-246. 

42.  Candes  E,  Wakin  M,  Boyd  S.  Enhancing  sparsity  by  reweighted  11 
minimization.  J.  Fourier  Anal.  Appl.  2008;  14:  877-905. 

43.  Stern  A,  Donoho  D,  Hoch  J.  NMR  data  processing  using  iterative 
thresholding  and  minimum  II -norm  reconstruction.  J  Magn  Reson 
2007;  188:  295-300. 

44.  Chen  SS,  Donoho  DL,  Saunders  MA.  Atomic  decomposition  by  basis 
pursuit.  Siam  J.  Sci.  Comput.  1998;  20(1):  33-61. 

45.  Kim  SJ,  Koh  K,  Lustig  M,  Boyd  S,  Gorinevsky  D.  An  interior-point 
method  for  large-scale  11-regularized  least  squares.  IEEE  J.  Sel.  Topics 
Signal  Process.  2007;  1(4):  606-617. 

46.  Barna  J,  Laue  E,  Mayger  M,  Skilling  J,  Worrall  S.  Exponential  sampling, 
an  aternative  method  for  sampling  in  two-dimensional  NMR 
experiments.  J  Magn  Reson  1987;  73:  69-77. 

47.  Eddy  M,  Ruben  D,  Griffin  R,  Herzfeld  J.  Deterministic  schedules  for 
robust  and  reproducible  non-uniform  sampling  in  multidimensional 
NMR.  J  Magn  Reson  2012;  214:  296-301. 

48.  Schmeider  P,  Stern  A,  Wagner  G,  Hoch  J.  Quantification  of 
maximum-entropy  spectrum  reconstructions.  J  Magn  Reson 
1997;  125:  332-339. 


NMR  Biomed.  2014;  27:  191-201 


Copyright  ©  2013  John  Wiley  &  Sons,  Ltd. 


wileyonlinelibrary.com/journal/nbm 


201 


SPLIT-BREGMAN-BASED  GROUP-SPARSE  RECONSTRUCTION  OF  MULTIDIMENSIONAL 

SPECTROSCOPIC  IMAGING  DATA 


Brian  Burns*,  Neil  Wilson \  M.  Albert  Thomas*' 

*  Department  of  Biomedical  Engineering,  UCLA,  Los  Angeles,  CA  90025  USA 
^Bio-Medical  Physics  IDP,  UCLA,  Los  Angeles,  CA  90025  USA 


ABSTRACT 

4D  Magnetic  Resonance  Spectroscopic  Imaging  data  pro¬ 
vides  valuable  biochemical  information  in  vivo,  however,  its 
acquisition  time  is  too  long  for  a  clinical  protocol.  In  this 
paper,  4D  simulated  MRSI  data  at  PSNRs=20,  10,  5  are 
under-sampled  2X,  4X,  6X,  or  8X  and  reconstructed  with 
Compressed  Sensing  and  Group  Sparsity.  A  Group  Sparse 
problem  solution  within  the  Split-Bregman  framework  is  de¬ 
rived,  which  allows  for  arbitrary,  over-lapping  grouping  struc¬ 
tures.  Preliminary  results  show  that  Group  Sparse  metabolite 
peak  amplitudes  and  line  width  reconstructions  are  more 
accurate  than  Compressed  Sensing  reconstructions. 

Index  Terms —  Group  Sparsity,  Compressed  Sensing, 
Split  Bregman,  Convex  Optimization,  Spectroscopic  Imaging 

1.  INTRODUCTION 

Magnetic  Resonance  Imaging  (MRI)  exploits  the  resonant 
frequency  of  1H  protons  within  water  in  vivo  to  generate 
anatomical  images  of  the  human  body.  Magnetic  Resonance 
Spectroscopic  Imaging  (MRSI)  is  a  similar  imaging  tech¬ 
nique  to  MRI,  however  in  lieu  of  the  resonant  frequency  of 
1  H  protons  in  water,  the  resonant  frequencies  of  1H  protons  in 
metabolites  are  used  to  generate  a  metabolic  image  instead. 
Each  metabolites  has  a  unique  resonance  spectrum  in  vivo 
caused  by  their  molecular  bonding  stmcture,  which  can  be 
used  to  identify  and  quantify  their  concentrations  within  each 
voxel.  Using  this  concentration  information,  the  biochemical 
compositions  of  healthy  and  diseased  tissues  can  be  deter¬ 
mined  without  the  need  for  invasive  biopsies  and  the  altered 
metabolism  of  cancers  can  be  detected  [1], 

The  Echo-Planar  Correlated  Spectroscopic  Imaging  (EP- 
COSI)  MRSI  pulse  sequence  allows  for  the  simultaneous 
acquisition  of  two  spatial  and  two  spectral  dimensions, 
(ky,kx,t-2,ti),  in  one  scan  in  vivo  [1],  This  4D  sequence 
provides  a  2D  spectrum  for  each  voxel  in  a  slice.  The  overlap 
of  resonance  peaks  within  a  single  spectral  dimension  is  a 
major  impediment  to  identifying  individual  metabolites  and 
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the  increased  spectral  dispersion  offered  by  a  second  spec¬ 
tral  dimension  can  disentangle  complex  over-lapping  spectral 
peaks.  However,  4D  MRSI  acquisitions  are  slow  compared 
to  most  MRI  sequences  and  can  take  up  to  40  minutes,  which 
is  too  long  to  be  used  on  a  routine  clinical  basis. 

4D  EP-COSI  data  acquisition  is  a  rasterized  scan  that  ac¬ 
quires  a  2D  spatiotemporal  plane,  kx—t-2,  from  the  4D  volume 
during  each  repetition  time  (TR).  The  second  spatial  (ky)  and 
spectral  (G)  dimensions  are  incrementally  acquired  between 
rasters  until  the  entire  4D  volume  is  sampled.  ky  is  incremen¬ 
tally  acquired  using  standard  phase  encoding  techniques  from 
MRI  and  G  is  incrementally  acquired  as  a  series  of  ID  spec¬ 
tra  with  modified  timings  [1],  To  accelerate  the  collection  of 
4D  EP-COSI  data  in  vivo  the  incrementally  acquired  dimen¬ 
sions,  ky  and  G>  can  be  non-uniformly  undersampled  (NUS). 
However,  NUS  produces  artifacts  that  must  be  removed  by  re¬ 
constructing  the  missing  samples  in  the  ky  — G  plane  through 
non-linear,  iterative  reconstruction  [2] . 

Previous  work  has  demonstrated  the  feasibility  of  under¬ 
sampling  the  mixed-domain  ky  —  t  \  plane  of  a  4D  MRSI 
data  set  and  reconstructing  the  missing  samples  with  Com¬ 
pressed  Sensing  (CS),  Total  Variation  (TV)  denoising,  and 
Maximum  Entropy  (MaxEnt)  [3,  4],  NUS  rates  as  low  as 
5X  were  demonstrated  in  vivo  and  showed  that  it  is  possi¬ 
ble  to  accelerate  the  acquisition  of  4D  MRSI  data  down  to 
a  clinically  acceptable  5-10  minutes.  The  current  work  uses 
the  Group  Sparse  (GS)  reconstruction  method  to  reconstruct 
simulated  4D  EP-COSI  data  sets  at  different  NUS  rates  and 
PSNRs  using  the  Split-Bregman  iterative  reconstruction  algo¬ 
rithm  [5,  6],  The  Split-Bregman  algorithm  has  been  modified 
to  allow  for  arbitrary  transform  coefficient  grouping  patterns 
of  overlapping  or  non-overlapping  groups  and  can  success¬ 
fully  reconstruct  NUS  MRSI  data  sets  down  to  8X.  The  Split- 
Bregman  algorithm  has  been  previously  used  for  TV-based  re¬ 
construction  of  NUS  MRSI  data  [3]  and  for  multi-channel  re¬ 
construction  of  NUS  MRI  data  by  extending  the  algorithm  to 
accomodate  row-wise  grouping  of  jointly  sparse  samples  [7], 

2.  THEORY 

2.1.  Split  Bregman  Algorithm 

The  Split-Bregman  algorithm  is  from  the  class  of  Alternating 
Direction  Method  of  Multipliers  (ADMM)  that  split  a  con- 


strained  problem  into  a  sequence  of  simpler  unconstrained 
sub-problems  [6],  It  differs  from  continuation-based  meth¬ 
ods  by  keeping  the  values  of  any  Lagrange  multipliers  fixed 
between  iterations  and  modifies  the  data  instead.  This  has 
the  benfits  of  increased  numerical  stablity  and  a  lower  depen¬ 
dence  on  the  initial  Lagrange  multiplier  values.  If  we  wish  to 
use  the  Split  Bregman  algorithm  to  solve  the  problem: 

min u,z  E(u,  z)  (i) 

such  that  H(u,z) 

where  E(u,  z)  is  a  convex,  non-differentiable  function  and 
H(u,  z)  can  be  assumed  to  be  of  the  form  Au  =  z.  We  first 
convert  it  to  an  unconstrained  problem: 

mmE(u,  z)  +  XH(u,  z)  (2) 

U,Z 

where  conventional  continuation-based  methods  would  in¬ 
crease  A  — >  oo  to  find  a  solution;  however,  we  apply  the 
Bregman  distance  relaxation  to  equation  2  and  split  it  into 
sub-problems  instead.  The  Bregman  distance  for  the  function 
E(u)  at  the  point  uk  is: 

DpE(u,uk)  =  E(u)  -  E[uk)  -  (pk,u-  uk)  (3) 

where  pk  is  the  subgradient  of  E(u)  at  uk.  The  Bregman 
distance  and  iteration  scheme  are  then  applied  to  equation  2: 

(uk+1,  zz+1)  =  min„>2  Dpe(u,  uu,  z,  zk)  +  A H(u,  z) 

pk+1  =  pku-X\7uH{uk+1,zk+1) 

Pkz+1  =  P\  -  A VzH{uk+\zk+l)  (4) 

which  is  iterated  over  k  until  convergence.  A  is  never  in¬ 
creased  and  under  fairly  weak  assumptions  on  E(u,  z )  and 
H(u,  z),  VUH  and  S7ZH  — >  0  as  k  — >  oo  so  that  pk+1  ->  pk 
and  the  Bregman  parameters  converge  [6].  However,  because 
H(u,  z)  is  of  the  form  Au  =  z,  equation  4  can  be  simplified: 

(uk+1,zz+1)  =  minM)Z  E(u,  z)  +  A||z  —  Au  —  bk ||| 

bk+ 1  =  bk  +  Auk+1  -  zk+1  ^ 

where  bk+1  is  a  Bregman  parameter  that  ensures  Au  —>  z  as 
the  iterations  converge  without  increasing  A  and  sacrificing 
stability.  Because  E(u,z)  is  convex  and  non-differentiable, 
equation  5  is  split  into  u  and  z  subproblems,  which  are  solved 
independently  at  each  iteration,  decoupling  u  from  z: 

uk+1  =  minu  E(u,  zk )  +  A||2fc  —  Au  —  bk ||| 

22+1  =  min,  E(uk  +  1 ,  z)  +  X\\z  -  Auk+1  -  bk ||| 

bk+1  =  bk  +  Auk+1  -  zk+1 

(6) 

2.2.  Group  Sparse  Reconstruction 

GS  reconstruction  is  an  extension  of  CS  that  exploits  the 
correlations  among  adjacent  transform  coefficients  caused 
by  their  structured  sparsity  [5],  Structured  sparsity  is  the 
tendency  of  large  transform  coefficients  to  be  adjacent  to 
eachother  and  form  clusters.  In  GS  reconstruction,  adjacent 
transform  coefficients  are  reconstructed  together  in  groups 


rather  than  individually,  as  is  done  in  CS.  By  reconstructing 
groups  of  coefficients,  the  GS  signal  model  correlates  indi¬ 
vidual  transform  coefficients  with  their  neighbors  allowing 
them  to  influence  eachother. 

GS  reconstruction  can  be  formed  as  a  constrained  convex 
optimization  problem  that  uses  the  /12-norm  as  the  objective 
function  instead  of  the  l  \ -norm  used  in  CS.  The  /12-norm  is: 

IN|l,2  =  ||zfll||2  +  ||zfl2||2  +  •••  +  ||^l||2  (7) 

where  zgi  €  Cp  =  { Uj,ui...un }  and  j,k,n  £  {5  = 

1...JV,  P  <  N}  is  a  group  of  transform  coefficients  from 
u.  The  MRSI  GS  reconstruction  problem  is  then  defined  as: 


min„jS  ||z||i,2 

such  that  || REu  —  /||  <  <r  (^) 

z  =  Gu 

where  u  =  (Y,  X.  F-2-  F\ )  is  the  reconstructed  spatial, 

spectral-domain  data,  T  is  the  4D  Fourier  operator,  R  is 
the  under-sampling  mask  that  determines  which  samples 
were  acquired  in  the  ky  —  t,\  plane,  /  =  ( ky ,  kx,t 2,  G)  is  the 
sampled  k-space,  time-domain  data,  G  is  the  group  matrix  of 
Is  and  Os  that  determines  which  coefficients  from  u  belong 
to  each  group  in  z  [8],  and  a  is  the  standard  deviation  of  the 
noise  in  /.  The  SB  problem  reduces  to  the  CS  problem  when 
each  group  contains  one  coefficient  and  G  =  I.  Because  each 
transform  coefficient  may  be  within  more  than  one  group  and 
separate  reconstructions  must  be  created  for  each  version  of 
that  coefficient,  the  set  of  transform  coefficient  groups,  z, 
may  contain  more  points  than  u. 

2.3.  Split-Bregman  Based  Group  Sparse  Reconstruction 


The  MRSI  GS  reconstruction  problem  can  be  solved  by  the 
Split-Bregman  algorithm  by  defining  equation  8  as  an  uncon¬ 
strained  problem: 

min  E(u,  z)  +  \H(u,  z)  (9) 


where 

E(u,z)  =  ||z||i,2  +p\\REu-  f\\l 
H(u,  z)  =  \\z-Gu\\l 


GO) 


By  following  the  process  defined  in  equations  4-6,  we  can 
derive  u  and  z  sub-problems  and  a  Bregman  parameter  update 
that  solves  the  unconstrained  problem  in  9: 

=  min„  ^11  REu  -  f  |||  +  A||^fc  -  Gu  -  bk  f2 
zz+ 1  =mmz\\z\\h2  +  X\\z-Guk+1 -bk\\l 
bk+1  =  bk  +  Guk+X  -  zk+1 

(11) 

where  /<  and  A  are  Lagrange  multipliers.  The  zk+1  sub¬ 
problem  is  non-differentiable,  however,  its  equivalent  prob¬ 
lem  can  be  solved: 


zk+1  =  min z  E  [IMI2  +  X\\zgi  -  (Guk+1)gi  -  (bk)gi \\%\ 
=  gshrink{Guk+1  +bk,j,G) 

(12) 


where  gshrink  is  group-wise  shrinkage  over  each  group  [8]. 
Because  the  uk+1  sub-problem  is  differentiable,  optimality 
conditions  can  be  derived  for  u  and  simplified  to: 

(liPR'RF  +  A G'G)uk+1  =  gF'R'p  +  A G'  ( zk  -  bk) 

(13) 

where  G'G  is  a  diagonal  matrix  with  each  (G'G)  u  the  num¬ 
ber  of  groups  that  contain  transform  coefficient  uk+  >  [8],  If 
each  transform  coefficient  is  in  the  same  number  of  groups, 
so  G'G  is  a  multiple  of  /,  the  left-hand  side  of  equation  13  is 
circulant  and  can  be  inverted  by  the  Fourier  transform: 

uk+ 1  =  F'K~XT  (uPR'f  +  A G'  ( zk  -  bk ))  (14) 

where  K  =  ( //  R'  R  +  A  G'G)  is  diagonal.  If  each  transform 
coefficient  is  not  in  the  same  number  of  groups  and  G'G  is 
not  a  multiple  of  /,  uk+1  in  equation  13  can  be  solved  for  by 
the  Gauss-Seidel  method. 

Iterating  the  sub-problems  and  Bregman  update  over  k  in 
equation  11  solves  the  unconstrained  problem;  however,  to 
solve  the  constrained  problem  in  equation  8,  an  outer  iteration 
is  completed  over  i  that  updates  /  based  on  changes  to  uk+1: 

f+i  =  fi  +  ^  _  RFuk+^  (15) 

The  number  of  iterations  over  k  is  application  dependent,  but 
for  the  MRSI  GS  problem,  we  use  15  iterations  over  k  for 
each  iteration  over  i  and  iterate  over  i  until  the  normalized 
residual  error  between  uk+1  and  /*  is  less  than  10e-5. 

3.  METHODS 

CS  and  GS  reconstructions  were  quantitatively  compared 
using  a  simulated  4D  EP-COSI  data  set  that  contained  in 
vivo  concentrations  of  the  metabolites:  Aspartate,  Choline, 
Creatine,  GABA,  Glutamate,  Glutamine,  Glucose,  Lactate, 
Myo-Inositol,  and  N-acetyl-aspartate.  The  2D  spectrum  of 
each  metabolite  was  simulated  using  the  GAMMA  NMR 
libraries  with  the  following  parameters:  100  t\,  1024  t2, 
TR/TE=  1.5s/30ms,  and  spectral  bandwidths  of  1250Hz  and 
2000Hz  along  F\  and  F2,  respectively.  Each  2D  spectrum 
was  line-broadened  by  10  Hz  to  simulate  in  vivo  conditions 
then  apodized  by  a  sine-squared  filter  along  t-\  and  a  skewed 
sine-squared  filter  with  skew  parameter  0.5  along  t2.  They 
were  then  copied  into  the  central  6x6  region  of  a  16x16  spatial 
grid  to  simulate  a  spatial  distribution. 

Noise  was  added  to  the  simulated  noise  free  4D  EP-COSI 
data  set  to  model  PSNRs  of  20,  10,  and  5.  The  SNR  was 
varied  by  adding  univariate  Gaussian  noise  to  the  noise-free 
real  and  imaginary  channels  of  the  4D  EP-COSI  data.  It  was 
under-sampled  8X,  6X,  4X,  and  2X  and  then  separately  re¬ 
constructed  by  CS,  GS  with  non-overlapping  groups  (GSi), 
and  GS  with  over-lapping  groups  (GS2).  The  group  sizes  used 
in  (GSi)  and  (GS2)  were  (ti,t2)  =  (4, 16),  which  were  em¬ 
pirically  determined.  The  CS  problem  was  solved  by  mak¬ 
ing  the  group  size  one  coefficient,  ( ti,t2 )  =  (1, 1).  Because 
the  noise  was  random,  each  SNR  was  simulated  and  recon¬ 
structed  50  times  per  NUS  rate. 


4.  RESULTS  AND  DISCUSSION 


2  3  4  5  6  7  8 

NUS  Rate 


Fig.  1.  The  mean  linear  correlation  coefficient  between  the 
fully  sampled  and  reconstructed  4D  EP-COSI  data  sets  after 
8X,  6X,  4X,  2X  NUS  at  PSNRs  of  20,  10,  5 

The  mean  linear  correlation  coefficients  between  the  re¬ 
constructions  and  the  fully  sampled  simulated  4D  EP-COSI 
data  set  are  shown  in  figure  1.  Values  were  calculated  from 
the  50  reconstructions  performed  at  each  PSNR  and  NUS  rate. 
Error  bars  are  not  shown  because  the  standard  deviations  were 
three  to  four  orders  of  magnitude  smaller  than  the  mean  cor¬ 
relation  coefficients.  The  linear  correlation  coefficient  is  a 
measure  of  the  strength  of  the  linear  relationship  between  two 
vectors;  two  vectors  are  strongly  correlated  for  values  close 
to  1  and  —1.  As  can  be  seen  in  figure  1,  both  the  GS  recon¬ 
structions  have  higher  mean  correlation  coefficients  than  the 
CS  reconstructions  for  all  PSNR  values  and  NUS  rates  tested. 
Additionally,  the  GS2  reconstructions  have  higher  mean  cor¬ 
relation  values  than  the  GSi  reconstructions  at  higher  NUS 
rates  and  this  difference  increases  with  PSNR.  This  suggests 
the  metabolite  peak  line-shapes  and  amplitudes  are  more  ac¬ 
curately  reconstructed  by  the  GS  reconstructions  than  CS. 

Figure  2  shows  a  representative  2D  spectrum  from  the 
simulated  4D  EP-COSI  data  set  at  PSNR=5.  The  left  shows 
the  fully  sampled  spectrum  followed  by  the  8X  NUS  spec¬ 
trum  in  the  next  column,  then  the  CS  and  GS  reconstruc¬ 
tions.  The  bottom  of  figure  2  shows  a  cross  section  of  each  2D 
spectrum  as  indicated  by  the  vertical  dotted  line  in  the  fully 
sampled  spectrum.  In  each  cross  section,  the  fully  sampled 
peaks  are  plotted  in  dark  gray  and  the  NUS  and  reconstructed 
peaks  are  in  black.  Both  GS  reconstructions  restored  each 
of  the  metabolite  peaks  to  their  fully  sampled  amplitudes  and 
widths,  as  well  as  reduced  the  noise  floor  below  that  of  the 
fully  sampled  spectrum.  The  CS  reconstructed  Glx  peak  is 
smaller  than  the  fully  sampled  peak  and  the  GABA  peak  am¬ 
plitude  is  larger,  while  the  NAA  peak  is  restored  to  its  origi¬ 
nal  amplitude  and  line-width.  This  indicates  a  higher  degree 
of  non-linearity  in  the  CS  reconstruction  than  the  GS  recon¬ 
structions,  which  may  have  adverse  effects  on  metabolite  con¬ 
centrations  estimates  that  rely  on  relative  peak  areas. 


Fig.  2.  (Top)  Representative  2D  spectrum  from  simulated  8X  NUS  4D  EP-COSI  data  set  with  a  PSNR=5  (Bottom)  ID  cross- 
section  indicated  by  dotted  line  on  fully-sampled  spectrum  showing  the  Glx,  GABA,  and  NAA  peaks  for  the  NUS  or  recon¬ 
structed  spectrum  in  black  and  the  fully  sampled  spectrum  in  dark  gray 


The  effects  of  the  reconstruction  on  the  GABA  cross  peak 
are  illustrated  in  table  1.  The  mean  amplitudes  and  Full- 
Width  Half  Maximums  (FWHM)  along  F\  for  the  metabo¬ 
lite  peak  with  respect  to  the  fully  sampled  values  are  given 
for  the  4X  and  8X  NUS  data  sets,  and  the  CS  and  GS  re¬ 
constructions  for  PSNRs=20,  10,  5.  Values  were  calculated 
from  the  50  reconstructions  as  in  figure  1.  As  was  seen  in 
figure  2,  the  CS  reconstructed  GABA  peak  is  larger  and  nar¬ 
rower  than  the  fully  sampled  peak  at  each  PSNR  and  NUS 
rate.  The  GS  reconstructions  were  better  able  to  reconstruct 
the  GABA  metabolite  peak  to  the  appropriate  amplitude  and 
FWHM.  The  GS2  reconstructions  show  a  trend  of  better  rel¬ 
ative  performace  to  the  other  reconstruction  methods  at  high 
NUS  rates  and  PSNRs,  which  was  also  seen  in  figure  1. 


PSNR=20 

PSNR=10 

PSNR=5 

4X 

8X 

4X 

8X 

4X 

8X 

NUS 

0.733 

0.612 

0.730 

0.608 

0.726 

0.601 

CS 

1.092 

1.043 

1.100 

1.035 

1.128 

1.035 

3 

GSi 

0.986 

0.944 

0.977 

0.933 

0.958 

0.894 

c 

gs2 

1.025 

1.01 

1.011 

0.987 

0.982 

0.936 

a 

cd 

NUS 

1.018 

1.128 

1.032 

1.121 

1.047 

1.112 

CS 

0.897 

0.913 

0.850 

0.891 

0.780 

0.842 

GSi 

0.958 

0.966 

0.980 

0.973 

1.008 

0.985 

a 

gs2 

0.973 

0.977 

0.975 

0.974 

0.975 

0.970 

§ 

Table  1.  GABA  cross  peak  amplitude  and  f)  FWHM  relative 
to  fully  sampled  values. 

The  Split-Bregman  GS  problem  converges  in  fewer  iter¬ 
ations  than  the  CS  problem,  but  takes  more  time  to  compute 
per  iteration  because  of  group-wise  shrink.  The  current  Mat- 
lab  implementation  of  GS  takes  2-3  times  longer  to  converge 
than  CS  and  is  a  current  topic  of  research. 

5.  CONCLUSION 

This  work  shows  that  GS  reconstruction  of  4D  MRSI  data 
sets  is  a  viable  alternative  to  CS-based  methods.  The  Split- 


Bregman-based  GS  algorithm  was  developed  and  evaluated 
on  simulated  data  sets  at  varying  PSNRs  and  NUS  rates.  GS 
reconstruction  results  demonstrated  better  metabolite  peak  re¬ 
production  and  lower  non-linearity  than  CS.  Further  work  is 
required  to  determine  the  optimal  grouping  strategy  under  dif¬ 
ferent  experimental  conditions. 
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